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Semi-analytical approach to criteria for ignition of excitation waves 
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We consider the problem of ignition of propagating waves in one-dimensional bistable or excitable 
systems by an instantaneous spatially extended stimulus. Earlier we proposed a method (Idris & 
Biktashev, PRL, vol 101, 2008, 244101) for analytical description of the threshold conditions based 
on an approximation of the (center-)stable manifold of a certain critical solution. Here we generalize 
this method to address a wider class of excitable systems, such as multicomponent reaction-diffusion 
systems and systems with non-self-adjoint linearized operators, including systems with moving crit¬ 
ical fronts and pulses. We also explore an extension of this method from a linear to a quadratic 
approximation of the (center-)stable manifold, resulting in some cases in a significant increase in 
accuracy. The applicability of the approach is demonstrated on five test problems ranging from 
archetypal examples such as the Zeldovich-Frank-Kamenetsky equation to near realistic examples 
such as the Beeler-Reuter model of cardiac excitation. While the method is analytical in nature, it 
is recognised that essential ingredients of the theory can be calculated explicitly only in exceptional 
cases, so we also describe methods suitable for calculating these ingredients numerically. 
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I. INTRODUCTION 

A. Motivation 


reaction-diffusion system, 


du 

~dt 


= D S + f(u), 


(1) 


An excitable system is a nonlinear active system that 
has a stable resting state , so that a weak, sub-threshold 
stimulus causes a straightforward return to the rest, but 
a stronger, over-threshold stimulus can produce a signif¬ 
icant, qualitatively different response. When such a sys¬ 
tem is spatially distributed, response to an over-threshold 
stimulus has the form of a propagating excitation wave in 
a shape of a non-decaying pulse, and one usually speaks 
about an excitable medium. A closely related concept 
is a trigger wave in a bistable medium : this takes place 
when the medium does not completely recover after a 
pulse but switches into a different steady state; trigger 
waves also often occur as idealizations of fronts or backs 
of excitation pulses. Excitable and bistable systems are 
widespread in nature and technology. Historically, the 
concept of excitability was first introduced in biology for 
nerve cells, and then applied to their electronic analogues. 
Later it was extended also to many other types of biolog¬ 
ical waves of signalling and in population dynamics, as 
well as such diverse physical situations as combustion and 
other chemical reaction waves, self-heating in metals and 
superconductors, phase transitions, domain wall move¬ 
ment in liquid crystals, nonlinear optics, surface boiling 
and laminar-turbulent transition in fluid flows, to name 
a few. See e.g. HHH3 for some literature on the topic 
It is often important not only to know that a particu¬ 
lar system can support a non-decaying propagating wave, 
but also to know what initial conditions can lead to it. 
The threshold character of the response of excitable and 
bistable systems, which characterizes already their local 
dynamics, gets much more complicated in the spatially 
extended context: the outcome of the localized pertur¬ 
bation will depend on its spatial and temporal charac¬ 
teristics as well as on its magnitude and modality. The 
conditions for initiation of propagating waves can be very 
important in practical applications. For instance, in the 
heart, excitation waves trigger coordinated contraction of 
the muscle and the failure of initiation can cause or con¬ 
tribute to serious or fatal medical conditions, or render 
inefficient the work of pacemakers or defibrillators la¬ 
in combustion, understanding of initiation is of critical 
importance for safety during the storage and transport of 
combustible materials [12]. In several key industrial pro¬ 
cesses, involving heat-generating elements, an important 
safety concern is the boiling crisis, or transition between 
a low-temperature and a high-temperature regimes |13j . 
which can proceed via trigger fronts m- 


B. Problem formulation 

We consider a formulation of the problems of initia¬ 
tion of propagating waves in terms of one-dimensional 


where u (x,t) :KxK-> is a d-component reagents 
field, d > 1, defined for x £ K. and t £ R+, vector- 
function f : —X describes the reaction rates and 

D £ R dxd is the matrix of diffusivity. We assume that 
this system has an asymptotically stable spatially uni¬ 
form equilibrium, called resting state, 

u(x,t)=u r , f(u r ) = 0, (2) 

and an orbitally stable family of propagating wave solu¬ 
tions of the form 

u(x, t) — u w ( X C w f S w )j 

u w (oo) = u r , u w ( oo) = u_, 

where u_ is also a stable spatially uniform equilibrium, 
f(u-) = 0, which may or may not coincide with u r . 
When u_ = u r the propagating wave solution is known 
as a propagating pulse, otherwise we shall call it a prop¬ 
agating front, and refer to u_ as the post-front equilib¬ 
rium; then u r may also be called the pre-front equilib¬ 
rium. In |3|, c w > 0 is a fixed constant, the wave prop¬ 
agation speed, and s w is an arbitrary constant, the pa¬ 
rameter of the family. Roughly speaking we assume that 
([3]) and ([ 2 ]) are the only attractors within the part of the 
phase space of (|T|) that is of practical interest, and we 
seek to find the boundary of the basins of attraction of 

U-w • 

In these terms, we seek to describe the localized (in 
space and time) perturbations of the resting state u,. 
which can lead to the propagating wave solutions u w . 
A localized perturbation will in fact typically generate 
two waves propagating away from the perturbed site in 
the opposite directions; this is obvious for perturbations 
that are even functions of x. With that in mind, we aim 
at classification of the solutions of the system ([I]) set on 
x £ [0,oo), t £ [0,oo), supplied with the following initial 
and boundary conditions, 

u {x, 0 ) = u 0 (a;) = u. r -f u s (x), x > 0 , 
Du x (0,t) = — I s (i), t > 0, (4) 

in terms of their behaviour as t — X 00: whether it will ap¬ 
proach the propagating wave solution (“ignition”) or the 
resting state (“failure”), as illustrated in fig. 0a,b)[56]. 
The functions u s (a;) and I s (t) are assumed to have a fi¬ 
nite support, u s (a;) = 0 for x > x s , and I s (t) = 0 for 
t > t s . 

A typical formulation is when only one of u s (-) and 
I s (-) is nonzero. If the dependence is on just one pa¬ 
rameter, then one speaks about threshold value(s) of the 
parameter, separating the two outcomes. When there are 
two parameters, one can talk about a threshold curve, or 
a critical curve , see fig. [lj(c) . The simplest and standard 
formulations are: 
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FIG. 1: (Color online) (a,b) Response to a below- and above-threshold initial perturbation in ZFK equation, given by formulas 
fll|4|5|6|64p . Parameter values: 8 = 0.13, I a = 0, x s = 2.10 for both, sub-threshold u s = 0.3304831 (a) and super-threshold 
u s = 0.3304833 (b) cases, numerics using central difference centered in space with step = 0.15 and forward Euler in time 
with step At = 0.01. Dash-dotted black lines: initial conditions, bold solid black lines: the critical nuclei, (c) The corresponding 
critical strength-extent curve, separating ignition initial conditions from decay initial conditions, (d) Sketch of a stable manifold 
of the critical solution for the ZFK equation. The critical nucleus is represented by the empty circle; the critical trajectories, 
constituting the center-stable manifold, are shown in thin solid black lines. The family of initial conditions is represented by 
the dasli-dotted lines. The bold solid black line is the critical trajectory with initial condition in that family. The sub-threshold 
trajectories are represented by the thin blue lines, while the thin red lines represent super-threshold trajectories. Note that the 
point where the initial condition intersect the center-stable manifold is shown as the filled circle. 


• “Stimulation by voltage”: 

l s (t) = 0, u s (x) = U s X(x). (5) 

That is, the initial condition is the resting state u r , 
displaced by the magnitude defined by the param¬ 
eter U s with a normalized spatial profile defined by 
X( x). In all specific examples we shall use simply 
a rectangular profile of a width a; s , 


for the critical curve (s S) Z7 S ); by analogy with the other 
case we shall call it the “strength-extent curve”. 

We find it convenient to formalize the initiation prob¬ 
lem as one posed on the whole real line i£l, 

<9u <9 2 u „. . , . , , , 

— = D-^- + f(u) + h(x,t), (x, t)etx K+, 

u(x, 0) = u r + u s (a;), h(:r,t) = 0 for t > t s . (9) 


X(:r) = H(a: s — x)e, (6) 

where H(-) is the Heaviside step function and e £ 
R d is a constant vector defining the modality of 
the perturbation. Then the critical curve can be 
considered in the plane (x s , U s ). 

• “Stimulation by current”: 

u s (a;) = 0, I s (t) = I s T (t). (7) 

That is, the initial condition is the unperturbed 
resting state, but there is a constant current in¬ 
jected through the left boundary of the interval, 
where I s defines the strength of the current and 
T(t) its normalized temporal profile. The most 
popular formulation is that of a rectangular pro¬ 
file of duration t s , 

T(i) = H(i a - t)e, (8) 

where the fixed vector e determines which 
reagent(s) are being injected. The critical curve 
will then be in the plane ( t Sl I s ). 

Historically, the “stimulation by current” formulation has 
been most popular with electrophysiologists and a stan¬ 
dard term for the critical curve (t s ,I B ) is the “strength- 
duration curve”. We are not aware of any standard term 


where the initial condition is an even continuation of the 
one in Q, 


u s (-x) = u s (x) 


U s 'X.(x), x > 0, 
t/ s X(— x), x < 0, 


( 10 ) 


and the boundary condition at x = 0 in Q is formally 
represented by the source term 


h(x,t) = 2I s T(t)S(x). (11) 


where <5(-) is the Dirac delta function. 


C. Aims 

Mathematically, the problem of determining the con¬ 
ditions of initiation of propagating waves in excitable 
or bistable media is spatially-distributed, nonstationary, 
nonlinear and has generally no helpful symmetries, so the 
accurate treatment is feasible only numerically. How¬ 
ever, the practical value of these conditions is so high 
that analytical answers, even if very approximate, are 
in high demand. Historically, there have been numer¬ 
ous attempts to obtain such answers, based on various 
phenomenological and heuristic approaches, e.g. SMS]- 
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The motivation for our present approach may be traced 
to the results by McKean and Moll |2Q] and Flores | 2Tj . 
who established that the boundary in the space of ini¬ 
tial data of 0 between the basins of attraction of 0 
and (|3| is a stable manifold of a certain “standing wave” 
solution, later also known as the critical nucleus. The 
critical nucleus is a solution of (jTJ) which is a bounded, 
non-constant function of x , independent of t and is un¬ 
stable, with one positive eigenvalue. The appearance of 
such a critical nucleus solution and its role is illustrated 
in fig. 0 if the initial data are very near the threshold 
between ignition and decay, the critical nucleus appears 
as a long transient before the outcome becomes apparent, 
and this does not depend on the sort of initial data, as 
long as they are near the threshold. This understanding 
has led to attempts to describe the critical conditions 
using Galerkin style approximations |22j . with analyti¬ 
cal answers obtainable by subsequent linearization |23| . 
This idea of the stable manifold of the critical solution 
has also been used to develop sophisticated numerical 
schemes for describing the critical conditions [24] . We 
have demonstrated that linearization of the stable man¬ 
ifold without any Galerkin projection but directly in the 
functional space produces surprisingly good reasults for 
a simple bistable model [25] . 

In the present paper, we seek to further explore and ex¬ 
tend the method of [25]. We focus on the case of stimula¬ 
tion by voltage and the strength-extent curve, leaving the 
stimulation by current and the strength-duration curve to 
other publications (note though that a simple case of the 
strength-duration curve was considered in [25]). In the 
case of stimulation by voltage, the mathematical problem 
is one about the basin of attraction of a dynamic system 
in a functional space. We investigate how the quality of 
approximation produced by our method depends on the 
parameters that define various test systems. Moreover, 
we investigate the feasibility of improving the accuracy 
by using a quadratic rather than a linear approximation 
of the critical manifold, and related problems. Finally, we 
extend the method to the case where there are no critical 
nucleus solutions. This is observed in multicomponent 
reaction-diffusion systems, where it has been previously 
demonstrated that instead of a critical nucleus, one has 
unstable propagating waves, such as critical pulses [28] 
or critcal fronts m 

The structure of the paper is as follows. In Section [TT] 
we describe the proposed analytical methods, including 
both the linear and the quadratic approximations of the 
critical manifold for the case of the critical nucleus, and 
the linear approximation for the case of moving criti¬ 
cal solutions, as well as the (rather standard) numeri¬ 
cal methods used in the study. Subsequent sections are 
dedicated to specific examples of application of the de¬ 
scribed method. In Section m we consider the one- 
component reaction-diffusion equation with cubic nonlin¬ 
earity, known as the Zeldovich-Frank-Kamenetsky (ZFK) 
equation, or the Nagumo equation, or the Schlogl model; 
this section recovers relevant results from [25] and fur¬ 


ther investigates the parametric dependencies and the 
quadratic approximation for a model of a propagating 
front. In Section [TV] the same is done to a piece-wise lin¬ 
ear analogue of the ZFK equation, known as the McKean 
equation. The piece-wise linearity of this equation means 
that some results can be obtained in closed form, where 
it was not possible in the ZFK case. Another special 
feature of this equation is that its right-hand side is dis¬ 
continuous, which presents certain technical challenges. 
The subsequent three sections are dedicated to exam¬ 
ples with moving critical solutions. Section [V] presents 
results in a two-component model where the critical so¬ 
lution is a moving front. It is a caricature model of 
cardiac excitation propagating fronts, and shares with 
the McKean equation the advantage of being piecewise- 
linear and of admitting analytical treatment, and also the 
challenge of having discontinous right-hand sides. Sec¬ 
tions [Vi] and [VH] are dedicated to two models where the 
critical solution is a pulse. Section [Vi] is an application 
to the FitzHugh-Nagumo (FHN) system which is a clas¬ 
sical “conceptual” model of excitable media, while Sec¬ 
tion [VlT] considers a variant of the detailed Beeler-Reuter 
(BR) ionic model of cardiac excitation, which, although 
not being physiologically precise from a modern view¬ 
point, includes many features of the up-to-date realistic 
cardiac excitation models. Both FHN and BR models 
do not admit full analytical treatment and we present 
the result of a hybrid approach, where the ingredients 
of the linearized theory are obtained numerically. We 
conclude by discussing the results and future directions 
in Section IVIIII 


II. METHODS 

A. Linear approximation of the center-stable 
manifold: principal approach 


We seek a classification of the outcomes in problem 
(9| depending on the parameters of the initial conditions 
(10), with h(a;, f) = 0. 


The principal assumption of our approach is the ex¬ 
istence of a critical solution , which is defined as a self¬ 
similar solution, 


u(s, t) = u(x — cf), 

f (u), 

u(oo) = u r , u(— 00 ) = u_ 


d 2 u du 

° = D dF + e d{ 


( 12 ) 


(where u_ may be different from u_ but in our examples 
u_ = u r when u_ = u,.) which, unlike the propagating 
wave u w defined by ([3]), is unstable with one unstable 
eigenvalue. Naturally, the speed c of the critical solution 
is also entirely different from the speed c w of the stable 
wave solution. 

Similar to the stable wave solution, there is then a 
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whole one-parametric family of critical solutions, 

u(x — ct — s), s G R. (13) 


Due to this translation invariance, this solution always 
has one zero eigenvalue. Hence its stable manifold has 
codimension two, whereas its center-stable manifold has 
codimension one and as such it can partition the phase 
space, i.e. it can serve as a boundary between the basins 
of different attractors. Our strategy is to approximate 
this center-stable space. In the first instance, we consider 
the following linear approximation. 

Let us rewrite the reaction-diffusion system (RDS) 0 
in a frame of reference moving with a constant speed c, 
so that u(x, t ) = u(£, t), £ = x — ct — s, r = t, 


du 

Ifr 


= D 


u(£,0) 


<9 2 u 

5C 2 


<9u 


f (u), 


= u r + u s (£ + s). 


We linearize this equation on the critical solution, 
which is stationary in the moving frame 


u(£,r) = u(£) + v(£,r). 
The linearization gives 

dv 9 2 v dv 

~ T> W + C ~di; + ^ v ’ 


v(£, 0) = u r + u s (£ + s) - u(0, 


where 


F (0 = 


df 

du 


u=fi(€) 


(14) 


(15) 


(16) 


is the Jacobian of the kinetic term, evaluated at the crit¬ 
ical solution. 

Equation (15) is a linear non-homogeneous equation, 
with a time-independent linear operator, 


5 T v = £v + h, £ 4 D _ +C _ +F(e) . (17) 

For the sake of simplicity, let us assume that the eigen¬ 
functions of C, 


£Vj(0 = A,v ; (0 


(18) 


are simple and form a basis in an appropriate functional 
space, and the same is true for the adjoint C + [57:. An¬ 
other assumption, which simplifies formulas and is true 
for all examples considered, is that all eigenvalues im¬ 
portant for the theory are real. We shall enumerate the 
eigenpairs in the decreasing order of A j , so by assumption 
we always have Ai > A 2 = 0 > A 3 > .... 

Then the general solution of problem (151 in that space 
can be written as a generalized Fourier series 


v (£’ t ) ='52 a i( T )' v j(€)- 

3 


The coefficients aj will then satisfy decoupled ODEs, 


day 

dr 


Xj a j 


( 20 ) 


where 


oj(0) = (w,(0 



the scalar product 


is defined as 


( 21 ) 



a T bd£, 


and W j are eigenfunctions of the adjoint operator, 

r) 2 B 

£+W, = A,W„ £+ = D t ^-c-+F t (C), (22) 
which are normalized so that 




The solution of (20) is 

Oj(t) = e XjT aj(0). 


(23) 


By assumption, Ai > 0, and due to the translational 
symmetry, A 2 = 0 , and the rest of the spectrum is as¬ 
sumed within the left half-plane. Hence the condition of 
criticality is 


ai(0) = 0. 


Using the definition of ai(0), we have, in terms of the 
original model, 



w(e + s)} = (w 1 (a 


u(£) - u r 


(24) 


This is an equation of the center-stable space, i.e. a tan¬ 
gent space to the center-stable manifold of the critical 
solution. Note that this tangent space is different for 
every choice of the critical solution as identified by the 
choice of s. 


B. Linear approximation of the strength-extent 
curve 

1. General setting 

Let us now consider the typical formulation, when the 
spatial profile of the initial perturbation is fixed and only 
its magnitude is varied, 


(19) 


u s (x) = U B X(x). 


(25) 
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Then ([24]) gives 

U s (Wi(0 | X(£ + a)) = (w x (01 u(£) - u r ) 


or 


C/s = (26) 

X>i(s) 

where the numerator J\f-[ is a constant, defined entirely 
by the properties of the model, 


A/”i = (Wx(^) | u(£) — u r ) 


(27) 


and the denominator T>\ depends on the shift s, 


p 1 ( s ) = (w 1 (o|x(c + s) 


(28) 


Hence to get the ultimate answer, we need an extra con¬ 
dition to fix the value of the shift s. 


2. The case of critical nucleus 



FIG. 2: (Color online) The sketch of a center-stable man¬ 
ifold of a moving critical solution. The critical solution is 
represented by the dashed bold black line; the critical trajec¬ 
tories, constituting the center-stable manifold, are shown in 
thin solid black lines. The family of initial conditions is rep¬ 
resented by the dash-dotted lines. The bold solid black line 
is the critical trajectory with initial condition in that family. 
The sub-threshold trajectories are represented by the thin 
blue lines, while the thin red lines represent super-threshold 
trajectories. Note that the point where the initial condition 
intersect the center-stable manifold is shown as the filled cir¬ 
cle. 


This is the case when c = 0, i.e. the critical solution is 
stationary, and moreover it is even in x. Then there is a 
natural choice of s = 0 prescribed by symmetry. It can 
also be motivated directly by considering the problem for 
x £ R as an even extension of the problem for x £ R+. 
In this case the position of the critical nucleus is fixed, 
there is no translational invariance, no associated zero 
eigenvalue, and we can consider the stable space, tangent 
to the stable manifold, as symbolised in fig. m , rather 
than center-stable manifold. 

That is, we have x = £, t = r, u = u, u(—£) = u(£), 
and (261 gives the explicit expression for the threshold 


OO 

0 


If, further, the stimulation is done by a rectangular per¬ 
turbation of the resting state, 


X(£)=H(x s -£)H(z s + £)e, (30) 


the symmetry £ i—> — £ and the previous “intuitively ob¬ 
vious” choice of s is not generally applicable. Recall that 
our approach is based on linearization, whereas the orig¬ 
inal problem does not, in fact, contain small parameters. 
In this formulation, the criticality condition depends on 
an “arbitrary” parameter s. We select an optimal value 
of the parameter, so as to minimize the error in the pre¬ 
diction. This is done based on a heuristic, motivated by 
the “skew-product” approach to the dynamics of systems 
with continuous symmetries, such as in I28| for symme¬ 
try with respect to shifts in R as in our present case, and 
also in [29, 30] for symmetry with respect to Euclidean 
motions in R 2 . This approach considers solutions u(x,t) 
of (jTJ) in the form 


U (x,t) = U (£,t), 

where £ = x — s{t), t = t, so u(-,r) describes the evo¬ 
lution of the shape of the wave profile in a frame of ref¬ 
erence which moves according to the law defined by the 
shift s(f), and the dynamics of the shift s(t) is determined 
from an extra condition, such as 


then we have 

OO _ 

/ W 1 (£) t (u(£) — u r ) d£ 

U s = ^ s -. (31) 

/W!(£)Ted£ 

o 

3. The case of moving critical solution 

This is the case when c > 0, and then we call the criti¬ 
cal solution either a critical pulse (for u_ = u r ) or a criti¬ 
cal front (u_ ^ u r ) 27], The problem now does not have 


= 0, (32) 

for an appropriately selected function /x(-), which allows 
to choose a unique value of s for any given profile u(-, t), 
from the class that is of interest to our study, at any given 
time t, perhaps with some inequalities to distinguish the 
front from the back. The specific examples of this extra 
condition, considered in [25H5D] included a condition on 
u(£,t) at a selected point £. A more generic way dis¬ 
cussed in [30] is to use any functional p(-) on u which is 
not invariant with respect to the group of translations of 
£, so that the functional could take a certain value only 
at selected values of s, say typically d s p, (u(£ + s, r)) ^ 0 
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and certainly d s fj, (u(£ + s ))| s=0 / 0. A popular and ef¬ 
ficient choice of such functional can be made when one 
considers perturbations of a relative equilibrium, as done 
e.g. in [5IH55| . This choice is based on the following ob¬ 
servation, adapted to our case of a one-parametric sym¬ 
metry group. An infinitesimal increment of the shift s is 
equivalent, in this situation, to a corresponding infinites¬ 
imal change of coefficient aj in an expansion like (19). 


Let Vj(£) = u'(£) be the “translational” mode, corre¬ 
sponding to A j = 0. Then a (locally) unique fixation of s 
can be achieved by requiring that aj =0. In our present 
situation, the index of the projector to the shift mode is 
j = 2. The resulting extra requirement is to be applied to 
the solution at all moments of time, including the initial 
condition, for which it gives 

(w 2 (£) u r + u s (£ + s)-u(£)^} = 0, 
leading to 

U s (w 2 (0 | X(£ + a)) = (w 2 (0 | u(£) - u r ) . (33) 

Another interpretation of the same requirement is that 
the condition a 2 (0) = 0, in addition to the already im¬ 
posed condition of criticality ai(0) = 0, makes sure that 
at least two first terms in the Fourier series (19) are zero, 
thus making v(£, 0) “smaller” in that sense. 

The two conditions give a system of two similar equa¬ 
tions, 


t/JMs) =M, 

U s T> 2 (s) = AT 2 , 


(34) 


where 


and 


= (w,(£) | u(£) -Ur), t 


= 1 , 2 , 


Z >/00 = (w,(£)|x(£ + s)), t 


= 1 , 2 . 


(35) 


(36) 


The definitions of M \, T >i here agree with the ones given 
earlier in (27), (28). We note that, if u_ ^ u r , inte¬ 
grals (35) converge if W^(£ -» — oo) —> 0 sufficiently 
quickly. 


System (341 is a nonlinear system of two equations for 
two unknowns, s and U s . It is linear and over-determined 
with respect to U s . The compatibility condition for the 
two equations for U s is N\D 2 (s) — Af 2 'D 1 (s) = 0, or 


*(0 X(C + S )) = 0 


where 


$(fl=MW 2 (0-J^W 1 (0, (37) 

presenting a nonlinear equation for s. For a rectangular 
stimulus profile, 

X(x) = H(x + x s ) H(a; s — x) e, 


the compatibility condition becomes 

— S+X s 

J *(0d£ = 0 

— s—x s 

where 

$(£) = e T <I>(£). (38) 

This equation for s can be transformed into a more con¬ 
venient form if we introduce the anti-derivative of 4>(£), 

Then 


V(s + x s ) - r](—s - x B ) = 0 , 


(39) 


that is, the two points £ + = — s + x s and = — s — x s 
are points of equal value of function ?y(-). If this func¬ 
tion happens to be unimodal, then a unique solution of 
the compatibility condition is guaranteed to exist, and if 
its monotonic pieces r] + {-) and are effectively in¬ 
vertible with, say, dom(? 7 + ) > pointwise, then 

(39) leads to a parametric equation for the critical curve 
U^x s ). If we denote the value of function rj(-) in (39) by 
C and take it as the parameter, then we have 


(40) 


e(o = (v ± r l (o, 
x s (o = \(t(o-r(o), 

s(o = -^(£ + (o+r(o), 

Ug {()=MM (8(0). 

For reference, we also summarise here the definitions of 


the ingredients of (40) given earlier: 


v(0=^iM0-^i(0, (4i) 

It(0 = J e r W e (?)d?, £ = 1,2 (42) 

OO 

Mt= J W[ (0 (u(0 - Ur) d£, £=1,2 (43) 


1 (D- 


(44) 


We note that in the case of a critical nucleus, c = 0, 
u is an even function, the operators C and C + commute 
with the operator of inversion £ i —> — £, the function Wj 
is even, the function W 2 is odd, AT 2 = 0, X 2 is even, i] 
is even, £ + = —, s = 0 and (40) formally recovers the 


result (261 obtained previously based on the choice s = 0 


as “intuitive” and “natural”. 


C. Quadratic approximation of the stable manifold 

The use of a linear approximation around the critical 
solution for the situation when distance from it is not 

























guaranteed to be very small is, admittedly, the weakest 
point of our approach. In this section, we consider the 
second-order approximation, in order to assess the limits 
of applicability of the linear approximation, and possibly 
to improve it. We restrict the consideration to the case of 
critical nucleus. We use the formulation oni€ (— 00 , 00 ) 
and on the space of even functions 

Rather than using the matrix notation as in the linear 
approximation, we shall now proceed with an explicit no¬ 
tation for the components of the reaction-diffusion sys¬ 
tems. We use Greek letters for superscripts to enumer¬ 
ate them, and adopt Einstein’s summation convention 
for those indices. In this way we start from the generic 
reaction-diffusion system 


du a 

Ik 


D af} 


d 2 uP 
dx 2 


fV), 


where the Fourier coefficients are defined by 

OO 

a j(t) = (Wj( x ) v(a;,t)^ = J Wj*(x) v a (x, t) Ax. 


Time-differentiation of this gives 

Oj = XjUj T ^ '] Qm,n^"ni a m 


(45) 


where 


y = O j 

&m,n y °6n,m 


Wf[x) f a 0 Ju(x))V£(x)V2(x) dx. (46) 


then consider the deviation v a of the solution u a from 
the critical nucleus u a , 

u a (x , t ) = u a (x ) + v a (x, t), 

the equation defining the critical nucleus, 

0 fpijP 

^+n*) = °, 

and the Taylor expansion of the equation for the devia¬ 
tion, 


We assume that eigenvalues are real and ordered from 
larger to smaller, Ai > 0 , A 2 = 0 is of course the eigen¬ 
value corresponding to the translational symmetry and 
an odd eigenfunction V 2 = u', and A j < 0 for all j > 3. 
Our task is to determine the conditions on the initial 
values of the Fourier coefficients 

OO 

Aj = dj( 0) = J Wj*(x)v a (x,0) dx 

— OO 

that would ensure that 


v a = uX + u)^^ + ..., 

where overdots denote differentiation with respect to 
time, subscripts (-) x denote differentiation with respect 
to space and Greek subscripts after a comma designate a 
partial differentiation by the corresponding reactive com¬ 
ponents. The right and left eigenfunctions are defined 
respectively by 

D af) d xx vf( x) + f$(x)Vf(x) = XjVfix) 


and 


Df* a d xx wf(x) + ff a (x)Wj (x) = A jWfix), 

where j G {1,2,3,...}, and the biorthogonality condition 
is 


OO 




Wf{ x)V?(x)dx = d j!k . 


— OO 


We consider only even solutions, so in subsequent sums 
only those j that correspond to even eigenfunctions are 
assumed. We seek solutions in the form of generalized 
Fourier series in terms of the right eigefunctions, 

v a {x,t) = Yl a i( t )Vj°‘(x) 

j 


ai(oo) = 0 , 

which means that the trajectory approaches the critical 
nucleus, so the initial condition is precisely at the thresh¬ 
old. 

Let us rewrite system (45) as an equivalent system of 
integral equations, 


a j (t)=e x i t (aj+ I e Xjt ' Q 3 m,n a m(t')a n {t') dt' 

\ q m,n 

Successive approximations to the solution can be ob¬ 
tained by direct iterations of this system, 

= 

Ai + / e~ Xjt ' J2 W } (f) dt' j . (47) 

q m,n J 

Taking a^ = 0 for all j, we have 
af ] = Aje Xjt . 

With account of Ai > 0, A j < 0, j > 3, and o^(t) —> 0, 
this implies that 




A\ = 0, 


Aj G K, j> 3, 
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which is the answer we have from the linear approxima¬ 
tion. The next iteration produces 


\t) = e> 


E 


Qrn,n.AmAr) 


3 ' / j \ _ \ _ \ 

'U / 'm S'n 

Tn,n>o J 


E Qm,nAmA n 

\ _ \ _ \ e 

/x 7 / 'm / 'n 

m,n >3 J 


Assuming that the sums converge, the last term always 
tends to zero as t —> oo because \ n < A 3 < 0 for all 
n > 3, and the first term tends to zero for all j > 3 for 
the same reason. So, the condition a^\t) —> 0 implies 
that that the first term vanishes for j = 1 , that is, 


A 1 


E 

m,n>3 


Qm,n-Am-A n 

Ai A m A n 


(48) 


which is our second-order (quadratic) approximation for 
the critical condition, as opposed to the first-order (lin¬ 
ear) approximation which states simply that A\ = 0. We 
see that the linear approximation will be more accurate 
when A n , n > 3 are smaller, and that for given magni¬ 
tudes of A n , the linear approximation will be better if 
Ai — A m — A n , the smallest of which is Ai — 2A3, is larger 
(remember we exclude all eigenpairs with odd eigenfunc¬ 
tions, including n = 2 ). 


Further iterations of (471 lead to still higher-order ap¬ 
proximations of the stable manifold of the critical nu¬ 
cleus, and possibly further improvement of the critical 
condition. This, however, is beyond the scope of this 
paper. 


Substitution into (48) of the definition of Aj in terms 
of the stimulation amplitude, 


OO 

Aj = J Wj(i) T (u r — u(a;) + f/ s X(a;)) da;, 


leads to a quadratic equation for U 8 , 

AU S 2 + BU S + C = 0 (49) 


where 


-A — ^ ^ Rrn,n^rn^n’> 

n,m>3 

B = — 2 'y ' Rm.n-AfmTAni (50) 

n,m> 3 

C = Ml ~\~ ^ ^ Rm,nMmMni 

n,m>3 


and 


Rm.n — 


Q 


i 

m.n 


Al A m A n 


— Rn,m i 


Afj = J Wj(x) T (u — Ur) da;, 

— OO 
OO 

Vj= J W 3 {x) T X(x)dx. 


(51) 


Note that the definitions of Afj, T>j here are the same as 
in (35), (36), with account of s = 0. 


An essential detail is the question of the properties of 
the spectra of C and C + . In the above derivation we as¬ 
sumed that these two spectra coincide, are discrete and 
all eigenvalues are simple. In the specific cases we con¬ 
sider below, these assumptions will be tested numerically; 
in particular, we shall observe that the spectra can in fact 
be continuous, so the formulas should be generalized, to 
replace summation over eigenvalues by integrals with re¬ 
spect to the spectral measure, and the convergence issue 
becomes even more complicated. 


D. Hybrid approach: numerical computation of 
functions required by the analytical theory 


1. Rationale 


The key to our linear approximation is the knowledge 
of u(a;), W^i) and, for thenon-self-adjoint cases, also 
of W 2 (t). For the quadratic approximation, ideally the 
whole spectrum of A^, W/-, Vg is needed. With a few 
fortunate exceptions, some of which will be discussed be¬ 
low, one does not have these analytically, so in practically 
interesting cases, one would need to employ a hybrid ap¬ 
proach, where these key ingredients are determined nu¬ 
merically before the analytical expressions (31) or (40) 
can be applied. The numerical problems can be posed as 
rather standard boundary-value problems, respectively 
nonlinear, for u, and linear, for A^, Wf, Vf. Here we 
describe the methods we used in specific examples pre¬ 
sented later. 


In all cases, for direct numerical simulation of time- 
dependent problems, we discretize the problems on a 
regular space grid on a finite interval x £ [0, L] as an 
approximation of 2 £ [0, oo), with fixed space step A x 
and a regular time grid with time step A t . Except where 
stated otherwise, we use second-order central difference 
approximations in space, with Neumann boundary con¬ 
ditions at x = L, and explicit first-order forward Euler 
method in time. 
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2. The case of critical nucleus 


a. Shooting Finding u means solving a nonlinear 
boundary-value problem. Most of the advanced meth¬ 
ods require a good initial guess for the solution. We find 
this initial guess by a version of the shooting method. We 
solve a sequence of the “stimulation by voltage” initial- 
value problem (5j6l, with fixed stimulation extent x s and 
varying amplitude U s . The sequence is selected with the 
goal of approaching the threshold value for U s which we 
denote as U*. This is done using the bisection method. 
Starting from an upper estimate U s , known to be suffi¬ 
cient for ignition, and a lower estimate t/ s , known to fail 
to ignite, we proceed with the following algorithm: 


repeat 

U* := -(Us + Us) (the new trial value of U s is 

the average of the current upper and lower estimates 
of the threshold); 

Solve the initial value problem with U s = U s * and 
determine if ignition or failure; 
if ignition then 

U s := U s * (the trial value ofU s will become the 
new upper estimate for the threshold); 

else 

U s := U s * (the trial value ofU s will become the 
new lower estimate for the threshold); 

end if 

until (|C/ s — C/ s | < tolerance) 


In fact, to achieve the best result, we typically use zero 
tolerance, i.e. repeat the bisection loop as long as U s * 
remains distinct from both U s and U s given the machine 
epsilon. The final value of U s * is the approximation 
of the critical amplitude U* for the given x s , the best 
achievable one at a given discretization. More precisely, 
we consider the last values of U s and U s as equally likely 
approximations U s * of U*, as which of them happens 
to be equal to ( U s + C/ s ) /2 in computer arithmetics is 
determined only by their position in the grid of floating¬ 
point numbers in the given architecture, rather than their 
relative merits as approximations. 

The so found estimate of U*(x s ) is used to determine 
an estimate for u. The idea is based on the observation 
that for U s very close to the exact threshold U*(x s ), the 
trajectory approaches the saddle point u(;r) to within 
a small distance and remains in its vicinity for a long 
time, see fig. 0 d )- So, to find the best estimate for the 
threshold trajectory obtained for U s = U s *, we calculate 


s(t) = Ml 2 



II W || 2 ( x,t ) dx 


(52) 


along the trajectory, find t* = argmin(S'(t)) and take 
n*(x) = u(x, t*) as an estimate of the critical nucleus 
u(:c). The result can be immediately used for the next 
step or serve as an initial guess for a more advanced 
boundary-value solver if a higher accuracy is required. 


Note that the key assumption of our theory is that the 
threshold manifold is the center manifold of a unique crit¬ 
ical nucleus solution, hence the above procedure should 
produce (nearly) the same u(cc) from any choice of x s . 
We used the procedure for different values of x s , both to 
verify the validity of this key assumption, and to assess 
the accuracy of the found critical nucleus. 

b. Marching The so found approximation of the 
critical nucleus u*(x) is then used to calculate the prin¬ 
cipal eigenvalue Ai and the corresponding eigenfunction 
W 3 . Since Ai is by definition the eigenvalue with the 
largest real part, we should expect that the solution of 
the differential equation 


<9w A -r<9 2 w -r, . 

St =c w = D s^ + F (I)W> 


(53) 


for almost any initial conditions, should satisfy 

w(x,t) = Ce Xlt Wi(i) (1 + o (1)), t —» oo. (54) 


for some constant C. We therefore consider the graph 
of ln|w(0,r)|, determine the linear part in it, and fit 
that linear part to a straight line by least squares; the 
slope provides an estimate of Ai. We have also veri¬ 
fied that the profile w(x,t)/w(0,f) remains virtually un¬ 
changed during this linear part, and took the most re¬ 
cent profile as Wi(i). Operationally, this is practically 
equivalent to the (more usual) procedure of estimating 
the eigenvalue Aj for a time interval from t to t + St 

as 8t~ l In ^w(a;, t + St) w(a;,f)^ / (w(x,t) w(x , t)^ 
and then ensuring that this estimate converges as t —> oo. 
Again, thus obtained Ai and Wi(i) can be immediately 
used or serve as an initial approximation for a more so¬ 
phisticated eigenvalue problem solver if a better accuracy 
is required. 

This method is of course easily extended to calculate 
not just the principal eigenpair, but a number of eigen- 
pairs with largest eigenvalues as long as they ar e real. 
If £ = C + then one only needs to calculate (53) for a 
number of linearly independent initial conditions, and at 
each step, in addition to normalization, also perform the 
Gram-Schmidt procedure. As discussed above, normal¬ 
ization of the first of the linearly independent solutions 
gives approximations of Wj and Ai. The second linearly 
independent solution is used to obtain a solution orthogo¬ 
nal to W i , which provides an approximation for W 2 , and 
its normaliztion provides and approximation for A 2 from 
its normalization. Then the third linearly independent 
solution is used to obtain a solution orthogonal both to 
Wi and W 2 , which gives an approximation for W 3 , and 
the corresponding normalization gives A 3 , and so on. In 
the non-self-adjoint case, C £ + , the orthogonalization 
will of course require this procedure to be done both for 
C and C + hand in hand, to calculate V 2 as orthogonal 
to Wi and W 2 as orthogonal to Vj and so on. 
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3. The case of moving critical solution 

a. Co-moving frame of reference We use the idea 
of symmetry reduction employed in Section [II B 3| for the 
theory, now for numerical simulations to reduce the prob¬ 
lem of moving critical solution to the case of a stationary 
critical solution, even though with non-self-adjoint lin¬ 
earization. To this end, we consider the equation 


<9u <9 2 u 

ffl =D &5 +f(u) - 


(55) 


where the position of the front, s, is defined implicitly by 
/i(u(s,f)) = 0, 

and perhaps some extra inequalities to distinguish the 
front from the back. Then in the comoving frame of ref¬ 
erence f = x — s(t), t = f, we have an unknown function 
of time and space, 

u(£,t) = u (x,t), 

and an unknown function of time, s(f), the system of 
PDEs and a finite constraint 

du „<9 2 u ds <9u 

a7 = D a? + d7af +f(u) - 


M(u(0 ,t)) = 0. 


(56) 


A relative equilibrium, including the moving critical so¬ 
lution, in the system (55), corresponds to an equilibrium 
in (56), so it is possible to use the same techniques de¬ 


veloped for the case of a stationary critical solution, to 
the comoving system (561. 

b. Shooting To find the critical solutions, we solve 


initial value problems for (56) and adjust the initial con¬ 


ditions so as to get as close to the initiation threshold as 
possible given the machine error. Solutions of (56) were 


found by semi-implicit time-stepping with the simplest 
(Lie) operator splitting, with four substeps, namely 

1. updating u by an explicit first-order (forward Eu¬ 
ler) scheme, for the nonlinear kinetics term f(u); 

2. updating u by a semi-implicit (Crank-Nicholson) 
scheme in time, and central-difference in space, for 

d 2 u 

the diffusion term D — 

ds 

3. finding the convection speed — based on a “vir- 

dr 

tual” or “predictor” convection substep, that would 
update u by an explicit in time, two-point upwind 
scheme without smoothing; 

4. the actual updating of u by an implicit in time, 
3-point upwind scheme with smoothing (Beam- 

ds 

Warming, 131]) for the convection term ———, us- 

dr 

ds 

ing the value of — found in the previous substep, 
dr 


As in the case of critical nucleus, the critical solution 
is estimated as the slowest point of the trajectory, i.e. 
at the moment r = r# = argmin ||u r \\ L2 . This includes 


ds 


both u#(£) = u(£,r#) and c# = — (t#). 

dr 

c. Continuation For the examples with non¬ 
stationary critical solutions, the accuracy of the critical 
solution found by shooting was often insufficient and we 
also found it as a solution of a boundary-value problem 
(12), which, in dynamical systems terms, is a problem 


of finding homoclinic (if u = u r ) or heteroclinic (if 
u_ ^ u r ) trajectories in a one-parametric family of 
autonomous systems for u(£), with parameter c. In 
the examples presented in this paper, we looked for 
homoclinics and used a simple and popular continuation 
method of finding such orbits, as a large-period limits of 
periodic trajectories of the same system, that is, 


Up dup 


0 = D 

u p {^ + P) = up (£) 


f (up), 


(57) 


using the continuation software AUTO Jj35|- When 
the problem is well posed, this defines a curve in the 
(P, Cp,up(£)) space. In our examples, the two ends 
of this curve extend to the limit P —> oo; one of 
the ends defining the stable propagating pulse solution, 
(c w , u w (£)), and the other defining the critical pulse solu¬ 
tion, (c, u(£)), which is of interest to us. An initial guess 
for the continuation procedure could be obtained from 
the shooting procedure described above, which would 
give the initial guess at the unstable branch of the curve, 
or just by direct numerical simulations of ([Tj) with P- 
periodic boundary conditions in x , which would give an 
initial guess at the stable branch of the curve. However, 
we preferred to use an ad hoc method, which is very 
popular for excitable systems, by finding the periodic 
orbits via Hopf bifurcation in a one-parametric exten¬ 
sion of (57), with an extra parameter corresponding to a 


“stimulation current”, that is an additive constant in the 
equation representing the dynamics of the activator, say 
the transmembrane voltage. 

AUTO uses collocation to discretize the solutions, and 
subsequent steps in our approach, such as marching and 
Arnoldi iterations, use u discretized on a regular grid. 
To interpolate the solution obtained by AUTO to the 
regular grid, we use piecewise Hermite interpolation [56] , 

d. Marching Once the critical solution u#(£) and 
its speed have been found, determination of the right 
and left eigenfunctions is done by calculating solutions of 
the initial value problems 

A=£v * d^ + c*A+F( £) v. (58) 


dr 


df 2 df 


and 


^C, 4D ^-^ + f(f).. (58) 
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The leading eigenvalue Ai and the corresponding right 
eigenfunctions Vi and left eigenfunctions Wj are ob¬ 
tained in the limit r —> oo for almost any initial condi¬ 
tions in (58) and (59). The second eigenvalue A 2 and the 


corresponding eigenfunctions V 2 and W 2 are obtained 
as linearly independent solutions of the same equations, 
satisfying the constraints 


W 2 Vi) = (W, V 2 ) =0, 


(60) 


using a Gram-Schmidt orthogonalization process, 
adapted to our non-self-adjoint situation. 

e. Arnoldi iterations When computing the required 
eigenvalues and eigenfunctions with the required accu¬ 
racy took too much time by the marching method dis¬ 
cussed above, we used the standard implicitly restarted 
Arnoldi iterations, using the implementation described 
in E3- This was applied to the matrices representing 
the right-hand sides of the discretized versions of equa¬ 
tions (59) and (58) to find left and right eigenfunctions re¬ 


spectively. We requested finding eigenvalues with biggest 
real parts and used the default values of the tuning pa¬ 
rameters. 


E. A priori bound in the critical nucleus case 


Finally we comment on a simple a priori bound for the 
critical curve, which follows from considerations different 
from the analysis of the central-stable manifold of the 
critical solution, and therefore may provide useful extra 
information. It applies for the case of d = 1, when the 
critical solution is the critical nucleus defined by 


du 

dt 


d 2 u 

dx 2 


+ f(u), 


(61) 


with the assumptions that f(uj) = 0, j = 1,2,3, u± = 
u r < 112 < U 3 , f(u ) < 0 for u £ (ui, U 2 ) and f(u) > 0 for 
u £ [ 112 , 113 ). In these terms, successful initiation means 
that at large t solution u[x, t) is a trigger wave from u\ to 
U3, and the failure of initiation means that u[x,t) —> u,\ 
as t —» 00 uniformly in x. 

It follows from the results by Fife and McLeod B38], 
that any initial conditions such that u(x, 0) £ [ 1 * 2 , 1 ^ 3 ] 
for x £ (— 00 , xi) and u[x, 0) £ [u\,U 2 ] for x £ (x 2 ,oo) 
guarantee ignition, and for rectangular initial conditions 
(25) this means that even for the smallest excess of U s 


over U 2 — Mi, this initial condition will produce ignition, 
provided that x s is large enough, so we have 


U*(x s )\U*, x s -^oo, 


(62) 


where 


U* = u 2 — u\. (63) 

In the following sections we verify this general method¬ 
ology by applying it to five examples. 


III. ZELDOVICH - FRANK-KAMENETSKY 
EQUATION 

A. Model formulation 


Our first example is the one-component reaction- 
diffusion equation, first introduced by Zeldovich and 
Frank-Kamenetsky (ZFK) [3H] to describe propagation 
of flames; it is also known as “Nagumo equation” [JD] 
and “Schlogl model” m- 

d= 1, D = (1), u = [u) , 

f ( u ) = (/(«)), f[u) = u(u-0){l-u), (64) 

where we assume that 6 £ (0,1/2). 

The critical nucleus solution u = (u) for this equation 
can be found analytically PUSS] Ei 


30^2 

(1 + 9)V 2 + cosh(a;\/0)v / 2 - 50 + 20 2 ' 


(65) 


The other two components required for the definition of 
critical curves in the linear approximation are Ai and 
Wi = Vi = (Vi) which are solutions of 


^- + (—3zt 2 + 2(0 + 1 )u -9)V 1 = AjVi, 

Ai > 0, Vi(±oo)=0. (66) 


We have been unable to find solution of this eigenvalue 
problem analytically. We note, however, that u given by 
(65) is unimodal, therefore v!, which is the eigenfunction 
of C corresponding to A = 0, has one root, hence by 
Sturm’s oscillation theorem, u' = V 2 and A 2 = 0, and 
there is indeed exactly one simple eigenvalue Ai > 0 and 
the corresponding V\ solving (66) has no roots. 


B. The small-threshold limit and the “fully 
analytical” result 


In this subsection we extend the results of [25] in the 
parameter space and correct some typos found in the 
latter paper. For 9 -C 1, the critical nucleus (65) is O ( 9) 
uniformly in x, and is approximately 

u[x) ~ -—- -j=r- = ^ 0sech 2 (ccV / 0/2). (67) 

' i + cosh(2;V0) 2 y ' > y > 

In the same limit, the nonlinearity can be approximated 
by f[u) ~ u(u — 9). With these approximations, problem 
(66) has the solution 


Ai » ^0, Vi » sech 3 (aV0/2), (68) 


and (31) then gives an explicit expression for the 


strength-extent curve in the form 

rr* 9?r0 

U* 


8 [4 arctan (e Xa ) + 2 tanh [x s ) sech (x s ) — 7r] 


(69) 
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(a) 


(b) 


(c) 


FIG. 3: (Color online) Strength-extent curves for the ZFK model, for 0 = 0.05, 0.15, 0.25, 0.35, 0.45 (bottom to top), comparison 
of direct numerical simulations (symbols) with theoretical predictions (dashed lines), (a) for the exact analytical answers in 
the # <g) 1 limit, linear approximation; (b) for t he h ybrid method, using the numerically found ignition eigenpairs, linear 
approximation; (c) same, quadratic approximation|70|. Discretization parameters: = 0.03, At = 4Aa, 2 /9, L = 100. 


where x s = \x s \/~9■ This approximation remains above 
the a priori lower bound (63), U* = 9, for all x s . 

Comparison of this approximation with the direct nu¬ 
merical simulations is shown in fig. d a )- We see that 
whereas for small 8 the comparison is reasonable for a 
wide range of x s , it quickly deteriorates at larger values 
of 0 , which is of course to be expected as the analytical 
expressions used are only valid in the limit of small 9. 


C. Hybrid approach 


Fig. a illustrates the processes of the numerical com¬ 
putation of the critical nucleus (a) and the ignition mode 
(b) in the ZFK model using the “shooting” algorithm 
described in Sectio n |II D 2[ for a selected value of the 
parameter 9. In fig. |4|a), the minimum of S(t) at about 
10 —5 , achieved at about t# « 50, designates the maximal 
proximity of the solution u(x,t&) of the nonlinear prob¬ 
lem (jT|) to the critical nucleus u{x), and so the former can 
be taken as an approximation of the latter. In fig. |4](b), 
the solution of the linear problem ( |53[ ), after an initial 
transient, mostly expiring before t = 10 , grows exponen¬ 
tially. The increment of this exponential growth gives the 
ignition eigenvalue Ai, and the corresponding solution 
profile, w(x, f)/w(0, t), which remains almost unchanged 
after t = 10, gives the ignition mode V\{x) = W\{x). 

The results of these numerical procedures are shown 
in fig. [5] We can see that the shooting procedure pro¬ 
duces good approximation of the critical nucleus, which 
for this case is known exactly, for all 9. We also see that 
the accuracy of the approximation obtained for ^ « 1 , 
unsurprisingly, is not good for larger 9 , see fig. |5[a). 
The solution of the adjoint problem shown in fig. pfb,c) 
demonstrates a nontrivial behaviour qualitatively differ¬ 
ent from the 9 < 1 analytical formulas: the eigenvalue 
starts deviating noticeably from ( 68 ) already for 9 ~ 0.1, 


and as 9 continues to increase across approximately 0.3, a 
qualitative change happens: the ignition eigenvalue Ai((9) 
stops increasing and starts decreasing, and the ignition 


mode V\ stops shrinking and starts expanding, and later 
even loses the unimodal shape and becomes bimodal, 
note the 9 = 0.45 curve in panel (b). The latter prop¬ 
erty should of course be expected: in the 9 /* 1/2 limit, 
the critical nucleus takes the form of two opposite look¬ 
ing fronts separated by the distance a ln(l/2 — 9), and 
the ignition mode is correspondingly a superposition of 
two sub-modes, each corresponding to its corresponding 
front, with the ignition eigenvalue A \ 0. 

One more observation can be made in fig. [5jc) about 
the behaviour of A y , j > 1. We see that the main assump¬ 
tions of the theory are satisfied and all these eigenvalues 
are negative, and moreover, they become more negative 
for larger 9. Further, the distance |A 3 — A 5 1 grows with 
0 , while a distance |A 5 — A 71 remains approximately the 
same and relatively small. This suggests that A 3 is a 
point of discrete spectrum, while A 5 and A 7 in fact rep¬ 
resent already the continuous spectrum and appear as 
discrete eigenvalues only due to the finite length L of the 
computational interval. This observation is confirmed by 
further study of these eigenvalues and the corresponding 
eigenfunctions: at increasing values of L, the distance 
| A 5 — A 7 I decreases, V\ and V 3 appear well localized to¬ 
wards the left end of the interval x £ [0, L], whereas 
V 5 and V 7 are manifestly non-localized, i.e. vary signifi¬ 
cantly throughout x £ [ 0 , 1 /] (not shown). 


Comparison of the resulting hybrid numeric- 
asymptotic prediction with the direct numerical 
simulations is shown in fig. [3](b). We see that for each 
value of 9 , the reasonable correspondence is observed 
in some range of x s . The large deviations are observed 
whenever U* gets large, which is expectable since the 
theory involves a linear approximation, and for large 
U* all Aj are large. We also note that for U* < 1, the 
quality of the hybrid approximation is in fact better for 
larger 9. This is also fully expectable based on the crud¬ 
est prediction of the quadratic theory: indeed one can 
see from fig. j5jc) that the spectral gap Ai — 2 A 3 , which 
is related to the accuracy of the linear approximation, 
grows with 9 (recall the discussion after equation (48)). 
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FIG. 4: (Color online) Illustration of the nu¬ 
merical computation of the critical nucleus 
and ignition mode by “shooting” and “march¬ 
ing” in ZFK. (a) Typical functions S(t) at 
near-threshold initial conditions in 0- 0’ 
(J5|, |6|, ( |64| ). Parameters: 9 = 0.15, x s = 0.6, 
U * « 1.1676..., m-m < 10 —5 , L = 20, 
A x = 0.02, A* = 4A;c 2 /9 . (b) Growth of the 
numerical solution of (531 in semilogarithmic 
coordinates, and its linear fit, defining the nu¬ 
merical value of Ai « 0.1425. 


S In |w(0, t)| 




(a) 


(b) 


V-i A 



FIG. 5: (Color online) Numerical computation of the components of the hybrid approach in ZFK. (a) Critical nucleus solu tion s 
for a 9 from 0.05 (bottom) to 0.45 (top) with step 0.1: numerical found by shooting, u*{x)\ exact analytical given by (651; 
approximate analytical for 9 1 given by ( |67| ). (b) The ignition mode, for a selection of values of 9, found by time marching 

based on numerical critical nucleus, (c) First four eigenvalues, found by marching based on numerical nucleus as functions of 9. 


D. Quadratic theory 


The quadratic theory result given by (49) involves dou¬ 


ble infinite sums over the stable modes of the linearized 
problems, so a practical application of this result in its 
fullness is problematic. However, we note that apart from 
the generalized Fourier coefficients of the critical nucleus, 
stimulus profile and the nonlinearity, this expression also 
has denominators increasing with the stable mode in¬ 
dices, so one may expect that depending on the proper¬ 
ties of the spectrum, the terms in the series may quickly 
decay and one can get a sufficiently accurate result by 
retaining only a few principal terms. As discussed in 
the previous subsection, for the ZFK equation, the lin¬ 
earized problem has one discrete stable eigenvalue and 
the rest of the stable spectrum is continuous. If we re¬ 
tain in (49) only the leading term, corresponding to the 


discrete eigenvalue, n = m = 3, we get a closed expres¬ 
sion for the critical curve, 

u * s » (70) 

-~T>i + \JR\ + 47?3 j 3 2?3 (A/i2? 3 — 

2 R^Dl 


the coefficients in which are defined by (46) and (511. 

The resulting approximations of the critical curves are 
shown in fig. |3](c) (to gether with the a priori bound 
U* = 9 given by (63)). Comparing those with panel 
(b), we observe that whereas there is little difference for 
larger 6 (the linear approximation for those was already 
reasonably good), there is noticeable improvement for 
6 = 0.05 and 9 = 0.15, where the quadratic correction 
term in (701 is more significant due to the relatively small 
denominator (Ai — 2A3). 


IV. MCKEAN EQUATION 
A. Model formulation 

Our second example is a piece-wise linear version of 
the ZFK equation, considered first by McKean in j401 
and then also in I4S1 : 

d= 1, D = (1) , u = (u) , 

f ( u )=(/0))> f(u) =-u + R(u-a), (72) 


or by expanding the square root, 

_ A/i Qls - V1N3) 2 
s ~ Vi V\ (Ai - 2A 3 ) 


where we assume that a £ (0,1/2), and H(-) is the Heav¬ 
iside step function. This model is a variant of ZFK, but 
with a special feature that makes it similar to the front 
model we consider later: the discontinuity of the kinetic 
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term. One of the practical issues caused by this disconti¬ 
nuity is that direct numerical simulations based on finite 
differences change qualitatively the behaviour of the sys¬ 
tem: the discretized critical nucleus solution, defined as 
an even, spatially nontrivial, stationary solution of the 
discretized equation, may not be unique and becomes 
stable, whereas in the differential equation it is unique 
and unstable. This phenomenon is akin to “propagation 
block” observed in discretized equations of the ZFK and 
McKean type and discussed e.g. by Keener [48] . with 
the exception is that here we are dealing with even so¬ 
lutions and spatially localized solutions (when extended 
to the whole line), as opposed to the trigger front solu¬ 
tions which are traditionally the object of interest in the 
context of propagation block. Keener’s result is about a 
generic system with smooth right-hand sides, and it pre¬ 
dicts “frozen solutions” for sufficiently large discretiza¬ 
tion steps. As we discuss in Appendix[Aj for the McKean 
model with its discontinous right-hand side, the situation 
is different in that the frozen solutions, at least formally, 
exist for all discretization steps, which motivates the use 
of a finite-element approximation, both in the direct nu¬ 
merical simulations for calculating the critical curve, and 
in the hybrid-method determination of the ignition mode. 
The finite element approach is discussed in Appendix |C| 
The critical nucleus solution in this equation is found 
exactly in a closed form, 

fl-(l-q) C0MX) 

u(x) = <^ 1 ; cosh(x*)’ 

[a exp(x» — x), 

where 

This solution is illustrated in fig. |6](a). 


x < x * 


(73) 


(74) 


B. Linear theory 


where 


2 a 


(^ e ” X * /a ) ’ 


^Wof-e- 1 */” 

zx* V a 


(77) 


and Wo(-) is the principal branch of the Lambert W- 
function as defined e.g. in [13] , The behaviour of 
this eigenpair at different values of a is illustrated 
in fig.[(|b,c). 

Substituting (731, (74), ([76]), ( [77] ) into (26), we obtain 


the analytical expression for the strength-extent curve, 


kA f 


x s < x* 


u:(x s ) = 


• 1 / \ ? *^s ^ u -'* 1 

smn(AvX s ) 


sinh(KX s ) — cosh(«;x s ) (e K ( x *~ x A — l) ’ 

x s > x*. 


(78) 


This prediction is compared with the direct numerical 
simulations in fig. [7](a). In this case, the theoretical pre¬ 
diction at larger x s falls below the apriori bound TJ* = a, 
so is “easily improved” by applying this bound. This is 
also shown in the figure. 

In this model, since the exact analytical solution for 
the critical nucleus and the ignition eigenpair is known 
for an arbitrary a £ (0,1/2), the “hybrid approach” is 
not necessary. For technical purposes we have tried it as 
well, and when used with finite-element discretization it 
works satisfactorily; but since it does not offer any extra 
insights, we do not present those results here. 


C. Quadratic theory 

The stable spectrum of the linearized problem in this 
case is entirely continuous, comprising all A < — 1 (see 
e.g. ED]), with the corresponding generalized eigenfunc¬ 
tions in the form 


Due to the discontinuity of the right-hand side, the lin¬ 
earization operator becomes singular, i.e. it contains the 
Dirac delta function (see Appendix [B] for a discussion): 

C = ^-2 - 1 — ~6(x - x*). (75) 

ox z a 

In the even-function extended problem, this is identical to 
the classical problem of a double-well potential in quan¬ 
tum mechanics. In the space of bounded functions, this 
operator has one positive eigenvalue. This eigenvalue and 
the corresponding eigenfunction can be written in the 
form 


Ai — —1 + K“, 

{ cosh (kx) 
cosh (kx*) ’ - 2 *’ 

exp (k(x* — x)), x > x*, 


(76) 


V(x; A) = p cos (px) — cos (px*) sin (p(x — x*)) H(x — x*) 


where p = \/—1 — A. 


Hence the sums in m,n in (49) are to be interpreted 


as integrals over A £ (—oo, —1]. The ensuing expressions 
are rather complicated and whereas it is plausible that 
the results can be expressed in a closed form, this goes 
well beyond the scope of this paper, and is left for an¬ 
other study. For now, as a proof-of-principle study, we 
have obtained a quadratic approximation of the critical 
curve, by restricting the infinite interval x £ [0, oo) to a 
finite interval x £ [0, L\, with a homogeneous Dirichlet 
boundary condition at x = L, thus making the spectrum 
discrete, and truncating the infinite sums in (49) to a 


finite number of terms. A represenative result is shown 
in fig. 0b). This was obtained for L = 10 and 287 eigen¬ 
values. 
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FIG. 6: (Color online) (a) Critical nucleus solutions, (b) ignition modes and 
various values of the parameter a. 


(c) eigenvalues of the McKean model (72 I for 


FIG. 7: (Color online) Strength-extent 

curves in McKean model: direct numerical 
simulations (red crosses) vs (a) linear the¬ 
ory, for a = 0.05 at the bottom increased by 
0.1 to a = 0.45 at the top, and (b) linear and 
quadratic theories, for a = 0.48. Blue long- 
dashed lines: analytical dependencies given 
by (781. Green short-dashed lines: the lower 
bound U s = a in (a) and the predictions 
given by quadratic theory in (b). Discretiza¬ 
tion: A x = 0.01, A t = 4 A;c 2 /9, L = 10. 
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V. THE CARICATURE MODEL OF THE 
/Na-DRIVEN CARDIAC EXCITATION FRONT 

A. Model formulation 

Our next example is the caricature model of an I^ a - 
driven cardiac excitation front suggested in [55] • It is a 
two-component reaction-diffusion system ([Tj) with u = 

(E, h) T , D = ^ ^ and f = (f E , fhf, where 

f E (E,h) = R(E~l)h, 

f h (E,h) = -(R(-E)-h), (79) 

T 

and H(-) is the Heaviside step function. The component 
E of the solution corresponds to the nondimensionalized 
transmembrane voltage, and the component h describes 
the inactivation gate of the fast sodium current, which is 
known in electrophysiology as Inh and which is mainly 
responsible for the propagation of excitation in cardiac 
muscle in the norm. 

A special feature of this model is that there is a con¬ 
tinuum of potential resting/pre-front states, 

u r = lim u = (—a, l) T , a > 0 , 

£->-o° 

and a continuum of potential post-front states, 
u_ = lim u=(w,0) T , uj > 1 , 

£—>■—00 


so any front solution connects a point from one contin¬ 
uum to a point from the other continuum. 

The critical solution u = (E, h) is described by 


m = 


uj — 


r 2 c 2 

1 + TC 2 


A/(tc) 


—a + ae 


-cf 


, £< A, 
£>-A, 


A/(tc) 


m = 


i, 


, £<o, 

£>o, 


(80) 


where the post-front voltage ui and the front thickness A 
are given by 


u = 1 + tc 2 (1 + a), A = - In 

c 


1 + a 


(81) 


and the front speed c is defined by an implicit equation 


c 2 In 


(1 + a)(l + tc 2 ) 


In 


a + 1 


= 0 , 


(82) 


or equivalently 


r = g(P,a) = ^/3- 1/<T , 


(83) 


where 


a = tc 2 , /3 


a/(a + 1 ). 


(84) 
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o’ 9/3 = o- go- = 0— 



FIG. 8: (Color online) Solutions of (831 for r = 7.7 and above 
with step 0.1. The dot is the global minimum of g(/3, a), the 
asterisks show the selected sets of parameters in the front 
model. 


Solutions of the transcendental equation (831 are illus¬ 
trated in fig.[ 8 j As shown in [42], f° r every r > t* ~ 7.674 
there is an interval of values of a, in which there are two 
solutions for c. The larger c corresponds to the stable, 
taller propagating front, and the smaller c corresponds to 
the unstable, lower propagating front. Our previous nu¬ 
merical simulations m indicated that the unstable front 
is the critical solution in this system: this unstable front 
is observed as a long-time transient for near-threshold 
initial conditions from different one-parametric families 
(corresponding to stimuli of different widths and varying 
magnitudes), which is a phenomenological evidence that 
its center-stable manifold has codimension one. 

We stress again that in this caricature model, the pre¬ 
front voltage —a is a parameter of the solution, rather 
than of the model, and for every r large enough, can 
take any value from an interval. In other words, each 
such resting state is not an isolated equilibrium, but is a 
member of a continuous, one-parametric family of equi¬ 
libria mo. 


The a priori bound (631 for the ignition threshold, 


discussed in Section |HE| is not applicable to the present 
two-component model. However, it is easy to see that if 
the initial condition E(x , 0) < 1 for all x, then the equa¬ 
tion for E reduces to a diffusion equation, and ignition 
of a propagating wave is therefore out of question. Since 
the initial perturbation is to be applied to the resting 
state E r = —a, we conclude that here there is a lower 
bound for the ignition threshold 


u: = i 


(85) 


which again can be used in conjunction with analytical 
predictions obtained from other considerations. 


B. Linearized problem and eigenfunctions 


The linearized operator © requires the Jacobian F(£) 
defined by (16). Formal differentiation of the Heaviside 
functions in the kinetic terms using the chain rule (see 


Appendix IB]) produces 
'dR(E-l) 


dE 


E=E 


3H(-A-£) /dE 

/ 


da 

1 


E'{- A) 


*(£ + A), 


dR(-E) 

dE 


E-E 


ma) 

da , 


'dE 

da 


E'{ 0) 


<5(0, 


then (79) gives for the Jacobian (16): 

(- 


F(0 = 


-<5(£ + A) H(-A-0\ 

E'(- A) 

1 


V 


-<5(0 ' 

tE'{ 0) T 




Hence the linearized equations (15) for v = ( E\,h \) are, 
componentwise, 


dE 1 

dr 


d 2 E 1 dE! 

da 2 +c ~da 


1 


E'{~ A) 

dhi 


dh\ 

ih =c ~da 


6 + A) hEi + H(—A - a) hi, 


— J 1 — S(a)E 1 --h 1 . ( 86 ) 

tE'{ 0) r 


The spatial operator in ( 86 ) is of the third order, so the 


linear eigenvalue problem (18) can be cast into a third- 
order ODE system for </>, % and ijj, where V = (</>,'!/’) 
and x = dcfi/da, which can be written in the matrix form 
as 




where E = (</>, x, V’) T an d 

( 0 1 

<5(£ + A) 


A = 


A + _ 

E'{- A) 

-m 

V tcE'( 0) 


0 \ 

-c —H(—^ — A) 

1 + At 


0 


/ 


(87) 


( 88 ) 


The regular part of matrix ( 88 ) is piecewise constant, 
hence the general solution to (87) can be written as 


H (0 = im exp (^m 0 > C e h 

m—1 


where symbol i takes one of three symbolic values, i = 
a , 6 , c, designating intervals I a = (oo, — A), If, = (—A, 0), 
I c = (0, oo), the vectors q^, q^, q c m , m =1,2,3 are 
the eigenvectors of A in these intervals, /i^, \i c m , are 
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the corresponding eigenvalues, and a m , b m , c m are coef¬ 
ficients of the solution in the bases of those eigenvectors 
in each of the intervals. We have, for A > 0. 


Mi = 


M 2 = 


1 + At 

TC 

—c — sjc 1 + 4A 


Ms = 


—c + \J c 2 + 4A 


= i'i(A) > 0, 

= -v 2 {\) < 0 , 
= ^3(A) > 0, 


for all three intervals i = a, &, c. For the sake of brevity, 
in the rest of this section we keep the dependence of ^ 1 , 2,3 
on A in mind, but omit in writing. 

Boundary conditions 3 (±00) = 0 [SO] then require that 
02 = ci = C 3 = 0, and by finding eigenvectors of A in 
the three intervals, we have 


X | = 


a 1 


h 


f 1 

Vl 

/ 0 > 
0 

V 1 ; 


pvii 


+ a 3 e 



^3? 





where 


q = + c) — A = ( ^ r N ) -I— >0. (90) 

TC J T 


These solutions are to be matched at £ = —A and £ = 0, 
with account of the singular terms in matrix ( 881 . Us¬ 


ing notation [• ] for a jump of a function at a point, the 
matching conditions for (87), ( 88 ), can be written as 


X(-A) 

V’(O) 


<K-A) 

E'{~ A)’ 

m 


tcE'{ 0) 

and we have continuity in all other cases, 

0(-A)l = [V’(-A)] = [ 0(0) 1 = [x(0) =0. 


For the solution (89), this amounts to the following alge¬ 


braic system for the coefficients: 
oie A + a 3 e "** A - b 2 e " 2 A - b 3 e ~ V3 A = 0 , 
ai acui e~ VlA + a 3 acv 3 e _ " 3A + b 2 e l ' 2A (acv 2 — e _i/A ) 
- b 3 e~ VsA ( acv 3 + e _i/A ) = 0 , 
ai q + bi = 0 , 
b 2 +b 3 - c 2 = 0, 
b 2 v 2 — b 3 u 3 — c 2 y 2 = 0, 

b 1 arc 2 + c 2 = 0. (91) 

where 

1 + re 2 


> 0 . 


TC 


(92) 


£ e l a , 


The solvability condition for this system is given by the 
roots of function f e ( A;...) (proportional to the Evans 
function) defined as 

f e (\;c,a,T) = ac(v 2 + v 3 ) e uA — 1 

,A -°' <*>) 

and the solution, up to normalization, is 



«1 

= 1 , 

«3 

= aTc 2 qe^ U2+V3,s>A - 

7 

(64 

bi 

= ~Q, 

0 


b 2 

= arc q , 

(6 4 

b 3 

= 0 , 


c 2 

= aTC 2 q. 

(89) 

The adjoint 

problem to (| 86 | is 


— e v 


(94) 


C+W = AW 


(95) 


where 


£ + =D T |^c^+F T (£), W =(<^) T , (96) 


and 


F t (£) = 


(97) 


-<5(£ + A) —<J(£)\ 

E'{- A) tE'( 0) 

1 

T 

Proceeding as in the previous case, we have a third-order 
system 


V H(-£-A) 


7 


d£ 

for 3 = (^>, x, V’) , X = <4, with the matrix 
/ 0 1 0 \ 


A = 


A + 3i+A) c -*(£) 


E'{- A) tcE'{ 0) 

-H(-£-A) 0 - (1 + Ar) 

C TC 


(98) 
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its piece-wise solution 

(l 


(m\ 

m 

\mJ 


= 


a 2 


bi 



h, 


-V3 | e-^^Gh, 
0 


ci 0 e 


-i'iS . 


c 3 


-^3 

0 


where 


P = 


c{v i + v 2 )' 


e~ V3i ,£, e / c , 

(99) 

( 100 ) 


the algebraic system for the coefficients stemming from 
the matching conditions, 


0,2 

e-^A _ 

- b 2 e-" aA - 

b 3 e^ A = 

0 , 

02 

acv 2 e 

-v 2 A 

- b 2 e 

v 2 A ( 

acv2 

+ e~ vA ) 


— b 3 

e^ A 

(-acy 3 + e 

— v A ^ 

= 0 , 

02 

pe~ V2 a 

'-fit 

e" 1 A = 

= 0 , 



b 2 

+ ^3 ~ 

C 3 = 

= 0 , 




b 2 

(XT CV 2 

- l>3 

a r cz /3 

+ ci 

+ c 3 

arch's = 

bi 

- Cl = 

0 , 






( 101 ) 


the solution of which, up to normalization, is 
a 2 = arc 2 (v 3 + v 2 ){v 2 + v 3 ) e (iyi+l,2) A , 

61 = arc (i / 2 + ^ 3 ), 

62 = — 1, 

b 3 = arc 2 (i/i + ^)(*z 2 + ^ 3 ) 6 ^-^) A + e-^+^l A , 

Ci = ar c (i/ 2 + ^ 3 ), 

c 3 = arc 2 (izi + v 2 ){u 2 + A + e~^ +V3) A - 1 , 

( 102 ) 

and its solvability condition is 


-(i/i+y 2 )A 

/ e = ac (i / 2 + v 3 ) H- 7 -r - e - " = 0 . 


- (103) 

TC{V 1 + I/ 2 ) 

The solvability condition |l03| is equivalent to ([93]), the 
characteristic equation of the linearized problem. Its de¬ 
pendence on A is implicit via v -\, v 2 , v 3 , p. The corre¬ 
sponding explicit expression is a available but too lengthy 
and we do not present it here. Alternatively, using sub¬ 
stitutions (84) and 


z = \J\ + 4A/c 2 , 


(104) 


Ve 



FIG. 9: (Color online) Plots of the charateristic function (1051 
for the four selected parameter sets. 



Set 1 

Set 2 

Set 3 

Set 4 

p 

0.5 

0.45 

0.4 

0.36 

T 

8.2 

8 

7.8 

7.7 

(7 

0.903152459 

1.036565915 

1.254739882 

1.559272934 

a 

1 

0.818181818 

0.666666667 

0.5625 

c 

0.331874289 

0.359959358 

0.401078655 

0.450003309 

Ai 

0.03990255 

0.035413196 

0.024380836 

0.006905681 


TABLE I: Selected sets of parameters and corresponding 
eigenvalues for the front model. 


the characteristic equation is explicitly rewritten in a 
more compact form, 


0 = v e (z) = ^<7 + ^1 + a ( Z ^ ^ 

+ (1 + z^) 2 / 4 - 1 . 



(105) 


Fig. [9 


illustrates the behaviour of the function v e (z) 


defined (105) for selected values of parameters /3 a nd a , 
indicated by asterisks in fig. § The roots z > 1 of ( |105[ ) 
define the positive eigenvalues Ai, and, of course, in all 
cases A 2 = 0 which corresponds to z = 1. Numerical 
values of Ai for the four selected sets of parameters are 
presented in Table |T] 

Knowing A^, £=1,2, we obtain the adjoint eigenfunc¬ 
tions Wf(£), £ = 1,2, by formulas (99). 


C. Strength-extent curve 


Given the expressions for the adjoint eigenfunctions 


(99), (102) for our model, we are now in a position to 
calculate the pre-compatibility function (41) required to 
obtain the analytical description of the critical curve. 
For the components of the left eigenfunctions, Wf = 
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FIG. 10: (Color online) Linearized theory ingredients for the four selected sets of parametes in the front model. Shown are 
components of vector functions, indicated in top right corner of each panel. The functions W, are scaled so that maximal value 
of /i-component is 1. Correspondence of lines with components is according to the legends at the top. 


(jt, ip e ^j , £ = 1 , 2 , we have 


for £ = 1,2. Further, (26) gives 


02 e " 2 S £e/ a , 

^(0 = ^e^ + ^e-^, i e h, 

£ e l c , 


V(s) = (W 1 (0 u s (£ + s )) = / ^(OdC, (106) 


c e 3 e 


and (42) gives 


( a? 2 f/ 

(€4, 

/• /• 

i>\0=\b[e-^, 

(€4, 

Z*(0 = / e T W,(C)d4= / 

[cfe-^, 

4 g 4c 

J J 


The resulting W^(£), £ = 1,2, for the selected values of 
parameters are shown in fig. | 10 [ against the correspond¬ 
ing critical nucleus solutions |61| . Then (351 gives 


a r t = (w,(o | u(0 - u r ) 


- e v )+- m d^ 


0 ^ 7 


ar 2 — 6 , _.A\ ,,e 


4e/ a 


^e-^ +V 3A + 

Vo Vo 


= 


= a 9 


+ bi 


uj + a — p e c(l + a)(p £ (l + rc 2 ) — r 2 c 2 )" 1 


1/9 tc + 1 


A 


a — (l + a)e " 2A a — (1 + a) e " 3 

— Uo 


i/fA 


_^_e-^S 

n c 




4 e J 6 , 


~r tr it it 
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3 g . i 
^ + c 


This general expression works for £ = 1; however for £ = 2 
it fails as A = 0 and consequently i/| = 0. The expression 
for £ = 2 can be obtained as the A —> 0 limit of the above, 
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or by redoing the integration for this special case. Either 
way, we get 


MO = 


~2 

tto j.2 t 

__p ^2 S 

o ^ 1 


+ 6 § (A + £), 


ALg—^2A _ ^|_ e -^ 2 A 

l/f ^2 

+6|A + dIC, 


a: 


I2 3 UTUJ 

Also, for A = A 2 = 0, Mi simplifies to 




u 2 c -„gA _ y2_ p -^A _f_ y2_ p i-|e 


g-^2^ + 


r e 2q 


Cel 6 , 


*2 


£el c . 


A4 = (w 2 (0|u(C)-u r ) 


_ 2 2arc 2 /_ 2 a ~ 2 \ . „ 

=a 2 “7T~—W + (a 2 - +&i'rc) 1 - 

c(l + rc^) V c / \ (1 + a)(l + tc 2 ) 

+ (b\ a + b\ \ A + 6 2 - + c 2 —. 

V / c c 

The critical curve is then described based on the func¬ 
tion r ](£), defined by ( |4l| ), 

7 7 (0=AA 1 T 2 (e)-A4Xi(C), 


using the implicit-function definition (40) (44). Fig. 11 


shows the function T](0 for the four selected parameter 
sets. It is clearly unimodal in all four cases, however the 
position of the maximum varies a lot. We note that for 
set 1 , the position of the maximum is far ahead of the 
critical front, which has a (rather unlikely) implication 
that this is where the main events, deciding whether the 
front will be ignited, take place. 

As can be seen from the above, even though the func¬ 
tion 7 /(£) for this example is found explicitly, its form ap¬ 
pears too complicated to establish analytically whether 
it is always unimodal, let alone to explicitly invert it. So 
we have done the inversion numerically. A comparison of 
the resulting critical curve against the direct numerical 
simulations is shown in fig. [12] We see that the accuracy 
of the theoretical predictions varies considerably, and for 
Set 4 is very close. We do not have any conclusive expla¬ 
nation of the difference in the accuracy, only note that 
better approximation is associated with a more reason¬ 
able position of the extremum of the pre-compatibility 
function ?)((;) and smaller values of Ai. 


VI. FITZHUGH-NAGUMO SYSTEM 
A. Model formulation 

The FitzHugh-Nagumo (FHN) system is a two- 
component reaction-diffusion system, which could be 


p 

C 

Ai 

^2 

0.05 

0.2561 

0.17204 

±1 

■ 10 -5 

0.13 

0.2328 

0.18619 

±1 

■ 10 -5 


TABLE II: Nonlinear and linear eigenvalues for the FitzHugh- 
Nagumo system. 


considered as a ZFK equation extended by adding a sec¬ 
ond, slow variable, describing inhibition of excitation. It 
is probably the single historically most important model 
describing excitable media. We consider it in the form 

d = 2, D = diag(l, 0), u = (u, v'j , 

f ( u ) = [fu(u,v),f v (u,v)^ , 

fu(u,v) = u(u - /3)(1 - it) - V, 

f v (u,v) ='y(au - v). (107) 


for fixed values of the slow dynamics parameters, 7 = 
0.01 and a = 0.37, and two values of the excitation 
threshold for the fast dynamics, /3 = 0.05 and /? = 0.13. 

Unlike the ZFK equation, the critical solution in sys¬ 
tem (107) is moving, as in the /N a front model, but it 
is a critical pulse rather than critical front. It is known 
(see e.g. [5B] and references therein) that in the limit 
7 \ 0 , this system has the critical pulse solution whose 
u-conrponent is small and u-component is close to the 
critical pulse of the corresponding ZFK equation. How¬ 
ever this does not provide good enough approximation 
for the linearized theory, and we used only the hybrid 
approach. We have obtained the critical pulse by nu¬ 
merical continuation of the periodic pulse problem using 
AUTO as discussed in Section [II D 3| the corresponding 
CV restitution curves are illustrated in fig. [13] For the 
critical pulses, we take the solutions at lower branches at 
P > 7.5 • 10 3 . The corresponding propagation speeds are 
given in Table [TT} 


B. Linear theory 


Fig. 14 and Table [XT] illustrate other ingredients re¬ 
quired for the semi-analytical prediction of the critical 
curves for the two selected cases. These are found by 
the straightforward marching method and then verified 
by Arnoldi iterations. We use discretization A x = 0.03, 
A t = 4 A,j, 2 /9, £ £ [— L, L], L = 100. As expected, |A 2 | 
are small. 

Fig. [I 5 ] shows the results of the calculation according 
to the formulas (40) (44). The “pre-compatibility” func¬ 


tion rj(£) defined by (41) in both cases is not unimodal, 
see fig. M a )> and, at least for /3 = 0.05, has two local 
maxima and one local minimum, hence implementation 
of the algorithm (40) (44) is not straightforward and re¬ 


quires investigation of the local extrema. We find that 
in both cases the adequate answer is given by the local 
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Set 1 Set 2 Set 3 Set 4 






FIG. 11: (Color online) Pre-compatibility functions r/(£) for the four selected sets of parameters in the front model. For 
visualization purposes, we show the normalized function, fj(£) = viO/viP) (dashed lines, right ordinate axes). For positioning, 
we also show the profile of the corresponding critical front (solid lines, left ordinate axes). 
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FIG. 12: (Color online) Strength-extent curves in the caricature 1 n 3 front model, for the four selected parameter sets. Shown 
also is the a priori lower bound for the ignition threshold given by equation (851. 
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FIG. 13: (Color online) CV restitution curves for the FHN 
model for two selected values of the model parameter. Stable 
(upper) and unstable (lower) branches are shown by different 
line types. 


maximum nearest to £ = 0, at the front of the critical 
pulse. The corresponding theoretical critical curves are 
shown in fig. EDM, in comparison with the curves ob¬ 
tained by direct numerical simulation. We observe that 
the theory works somewhat better for p = 0.05 than for 
P = 0.13, although the eigenfunctions shown in fig. [14] 
for the two cases look rather similar. Again, the bet¬ 
ter accuracy of the linearized theory here is associated 


with smaller value of Ai, although the relative difference 
between the two cases is small in itself. 


VII. MODIFIED BEELER-REUTER MODEL OF 
CARDIAC EXCITATION 

A. Model formulation 


Here we look at a variant of the classical Beeler- 
Reuter (BR) model of mammalian ventricular cardiac 
myocytes [46), modified to describe phenomenologically 
the dynamics of neonatal rat cells pTFHSO) : 


d= 7, D = diag(l, 0,0,0,0, 0, 0), 

u = (v,h, j,x 1 ,d,f,Gaj , 


f(u) 


^ —(iKi+Ixl+lNa + Is) ^ 

<*/>(! -h)~ p h h 

a x i(l - xi) - PxiX\ 
a d { 1 - d) ~ Pdd 
«/(!- /)- Pff 
\-10~ 7 I s + 0.07(10 -7 - Ca) j 


(108) 

(109) 


( 110 ) 
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FIG. 14: (Color online) FHN theory ingredients for (a) /3 = 
0.05 and (b) /3 = 0.13. Shown are components of scaled vector 
functions, indicated in top right corner of each panel, where 
u = Su, Wj = S~ 1 Wj, and S = diag(l, 10). The space 
coordinate is chosen so that £ = 0 at the maximum of u. 
Correspondence of lines with components is according to the 
legends at the top. 

The detailed description of the components of f(u) is 
given in Appendix [D] Important here is only the depen¬ 
dence on the “excitability” parameter ct, which appears 
in the equations in the following way: 

I Kl = 0.35 (0.3 - a) Tk~(V). 

In [SI], special attention was given to a = 0.105 (“less 
excitable”, with negative filament tension of the scroll 
waves) and a = 0.115 (“more excitable”, with positive 
filament tension of the scroll waves). These are also the 
two selected cases for our study here. 

As in the FitzHugh-Nagumo system, we obtain the 
CV restitution curves by continuation using AUTO, and 
use a solution at its lower branch as the critical pulse; 
see fig. [16] The corresponding propagation speeds are 
given in Table |III| 



-40 -20 0 20 40 -40 -20 0 20 40 f 

(a) (b) 

u; 



(c) (d) 

FIG. 15: (Color online) Results of linearized FitzHugh- 

Nagumo theory for (a,c) /? = 0.05 and (b,d) /? = 0.13. (a,b): 
The pre-compatibility function ry(£), used to compute the the¬ 
oretical critical curve. The u(£) component of the critical so¬ 
lution u is also shown for positioning purposes. (c,d): Com¬ 
parison of the theoretical critical curves obtained in the linear 
approximation, and the critical curves obtained by direct nu¬ 
merical simulations. 


c 



FIG. 16: (Color online) CV restitution curves for the BR 
model for two selected values of the model parameter. Stable 
(upper) and unstable (lower) branches are shown by different 
line types. 


marching method and then verified by Arnoldi itera¬ 
tions. We use discretization A x = 0.03, A t = 4 A;e 2 /9, 
£ £ [— L, L], L = 90. Again, we note that numerically 
found values of | A 2 1 are small. 

The pre-compatibility functions ry(£) (see fig. [I8|)a)) 
are this time nearly unimodal with a prominent maxi- 


13. Linear theory 


Fig. 

quired for the semi-analytical prediction of the critical 
curves for the two selected cases. As for the FitzHugh- 
Nagumo system, these are found by the straightforward 


17 and Table |III| illustrate other ingredients re- 


a 

c 

Ai 

A 2 

1.05 

0.04232 

0.01578 

00 

1 

0 

r—H 

<M 

-H 

1.15 

0.04497 

0.01515 

±1■ 10 -8 


TABLE III: Nonlinear and linear eigenvalues for the modified 
Beeler-Reuter model. 
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-60-40-20 0 20 40 60 -60-40-20 0 20 40 60 £ 

(a) (b) 

FIG. 17: (Color online) BR theory ingredients for (a) 

a = 0.105 and (b) a = 0.115. Shown are components 
of scaled vector functions, indicated in top right corner of 
each panel, where ... u = Su, Wj = 10 4 S“ 1 W j , and 
S = diag(10 -2 ,1,1,1,1,1,10 5 ). The space coordinate is cho¬ 
sen so that £ = 0 at the maximum of V. Correspondence of 
lines with components is according to the legends at the top. 


mum near the front or the peak of the critical nucleus. 
Again, despite apparent similarity of the eigenfunctions 
in fig. 17 between the two cases, the shape of the rj(£) 
graphs in fig. [l8|(a,b) is considerably different, and the 
resulting theoretical critical curves, shown in fig. [l8|(c,d) 
are much better for a = 0.105 than for a = 0.115. This 
time, Ai is nearly the same in both cases. 


VIII. DISCUSSION 

In this paper, we have substantially extended the 
method proposed in [25] for an analytical description of 
the threshold curves that separate the basins of attrac¬ 
tion of propagating wave solutions and of decaying so¬ 
lutions of certain reaction-diffusion models of spatially- 
extended excitable media. The method is extended in 
two ways. Firstly, it is generalized to address a wider 
class of excitable systems, such as: 

• multicomponent reaction-diffusion systems; 

• systems with non-self-adjoint linearized operators; 



FIG. 18: (Color online) Results of BR theory for (a,c) 

a = 0.105 and (b,d) a = 0.115. (a,b): The pre-compatibility 
function ??(£), used to compute the theoretical critical curve 
(dashed lines, right ordinate axes). The V(f) component of 
the critical solution u is also shown for positionining purposes 
(solid lines, left ordinate axes). (c,d): Comparison of the the¬ 
oretical critical curves obtained in the linear approximation, 
and the critical curves obtained by direct numerical simula¬ 
tions. 


• in particular, systems with moving critical solutions 
(critical fronts and critical pulses). 

Secondly, the method is extended from being a linear 
approximation to being 

• a quadratic approximation 

of the stable manifold of the critical nucleus solution, re¬ 
sulting in some cases in a significant increase in accuracy. 

The essential ingredients of the theory are the criti¬ 
cal solution itself, and the eigenfunctions of the corre¬ 
sponding linearized operator. For the linear approxima¬ 
tion in the critical nucleus case, we need the leading left 
(adjoint) eigenfunction; in the moving critical solution 
case, we need two leading left eigenfunctions; and for 
the quadratic approximations we need as many eigen¬ 
values and left and right eigenfunctions as possible to 
achieve better accuracy. Of course, closed analytical for¬ 
mulas for these ingredients can only be obtained in ex¬ 
ceptional cases, and in a more typical situation a “hy¬ 
brid” approach is required, where these ingredients are 
obtained numerically. Still, we believe that this approach 
offers advantage over the determination of the excitation 
threshold by direct numerical simulations, both in terms 
of insight and computational cost. 

It is still an open question, why in some cases our 
method works better than in others. A partial answer 
to that question is offered by the quadratic approxima¬ 
tion: the linear approximation performs better when the 
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corrections offered by the quadratic approximation are 
small. This seems to work for scalar equation with sta¬ 
tionary critical solutions (critical nucleus). However, in 
this paper we have not investigated the quadratic ap¬ 
proximation for the cases of moving critical solutions, 
that is, critical fronts and critical pulses. Here one has 
to bear in mind that the theory was presented here under 
the assumption that the spectrum is real, whereas in the 
non-self-adjoint case it does not have to be. So, this ques¬ 
tion remains an interesting direction for future research. 
Progress in that direction may help to understand when 
the linear approximation works better in such cases. 
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Appendices 


Another direction for future research is the possi¬ 
ble extension of the theory to different initiation pro¬ 
tocols, most notably, to the case of strength-duration 
curves for stimulation of an excitable cable by a stim¬ 
ulus localized in space and extended in time. In i23| we 
have shown that in the scalar case, the linearized theory 
readily gives the classical Lapique-Blair-Hill exponential 
rheobase-chronaxie expression. It is known that in realis¬ 
tic excitable systems, this formula does not always work 
well, and it is likely that the more complicated expres¬ 
sion coming out of our theory based on moving critical 
solutions would perform better. 

An obvious extension of our approach that is required 
for many applications is the extension to two and three 
spatial dimensions. 

It is also of interest to investigate whether the proposed 
semi-analytical approach to ignition of excitation waves 
can be adapted to address the reverse problem of estab¬ 
lishing conditions for decay (block) of an already prop¬ 
agating excitation wave. This question is of particular 
importance in practical situations, for instance for wild 
fire extinction, cardiac defibrillation and others. Some 
crude criteria for block of excitation can be established 
from asymptotic considerations of conditions when prop¬ 
agating wave solutions cease to exist, e.g. ESS, however 
extention of the approach presented in this paper may 
offer more refined criteria. Propagating waves do not 
have the shape of rectangular pulses, as typical stimuli 
do, and decay from a general wave form must be consid¬ 
ered in greater detail. An extra feature in two and three 
dimensions is the possibility of “wave breaks” which is a 
situation distinct from a complete decay and which is of 
particular relevance for cardiac arrhythmias. 


Appendix A: On “frozen nuclei” in the McKean 
equation 


Direct numerical simulations in this model are more 
difficult because in the standard finite-difference dis¬ 
cretization, the critical nucleus solution, defined as an 
even, spatially nontrivial, stationary solution of the dis¬ 
cretized equation, may not be unique and is stable, 
whereas in the differential equation it is unique and un¬ 
stable. This phenomenon is akin to “propagation block” 
or “propagation failure” observed in discretized equa¬ 
tion of the ZFK and McKean type and discussed e.g. 
in 03 EH El, with the exception that here we are deal¬ 
ing with even solutions, which correspond to spatially 
localized solutions when extended to the whole line, as 
opposed to the trigger front solutions which are tradi¬ 
tionally the object of interest in the context of propa¬ 
gation block. Keener’s 05] result is about a generic sys¬ 
tem with smooth right-hand sides, and it predicts “frozen 
solutions” for sufficiently large discretization steps. As 
we shall see, for the McKean model with its disconti- 
nous right hand side, the situation is different in that the 
frozen solutions, at least formally, exist for all discretiza¬ 
tion steps. 

For the McKean model, the “discrete critical nucleus” 
solutions can be studied analytically. For the regular grid 
discretization, 


Uj = u(xj), Xj = jA x , j £ Z, u-j = Uj, 


we have 


duy 
d t. 


Uj—i — 2 Uj T Uj .)_i 

A,. 2 


Uj + H(uj — a). 


Finally, we note that the problem of initiation of waves 
is of importance in all excitable systems, not just in car¬ 
diology. The theory presented here is likely to face new 
challenges in new applications. For instance, combus¬ 
tion waves sometimes can propagate in oscillatory man¬ 
ner, i.e. as relative periodic orbits [52], which makes it 
plausible that the critical solution there also is a rela¬ 
tive periodic orbit, and the transition to turbulence in 
shear flows, although exhibiting features of excitability, 
is in terms of models beyond reaction-diffusion even in 
the simplest phenomenological description m- 


We use notation A , B for the set of all integers j such 
that j > A and j < B, that is, A, B = [A, B] D Z. For 
the critical nucleus solution, we expect that there exists 
an integer m such that Uj > a for j £ — to, to and Uj < a 
otherwise. We ignore the exceptional case when Uj = a 
exactly for some j as it is not interesting in practice. 
Let Uj = Vj + H(uj — a). Then, separately on each of the 
intervals j £ — oo, — m — 1, j £ —to, to and j £ to + 1, oo, 
we have 

Vj -1 — 2 Vj + Vj +1 — A x 2 Vj = 0. (Al) 
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For Vj oc a 3 , this gives 

(T 2 — (2 + A x 2 ) <7 + 1 = 0, 
and so a = /i or a = 1//J,, where 


/r — 1 + A x /2 + A x y 1 + A x /4 > 1. 


(A2) 


We note that equation (A1) applied for j £ —to, m in fact 


involves Vj for j £ — m — l,m + 1, and the same equa¬ 
tion applied for j £ to + 1, 00 describes v 3 for j £ m, 00 , 
so the non-overlapping sub-intervals of the equation cre¬ 
ate overlapping sub-intervals in the piecewise described 
solutions. 


Considering (A1) for j £ m + 1 ,00 with account of the 


boundary condition lim^oo Uj = 0, we have 


Uj = Vj = Afi 


—j 


J £ TO, OO 


for some constant A. Further, considering (Al) for \j\ < 
to, we get the even solution in the form 


= 1 + Vj = 1 + B + [X J ) , j £ —to — 1, to + 1 


for some constant B. The matching condition to deter¬ 
mine A and B is that the two solutions should coincide 
at the overlap points, j = m and j = m + 1. This gives 

a A"*"; 1 - 1 ) _ *, 

fi m (fi + 1) fi m {n + 1 ) 


and the nontrivial time-independent solution Uj(t) = Uj 
in the form 


a 



FIG. 19: (Color online) Non-uniqueness of the discrete critical 
nucleus solutions is observed for some a > ai(A x ), and for all 
a > a 2 ( A x ). 


Considering that the discrete solution (A3) approximates 


the exact critical nucleus solution (731, we have the cor¬ 


responding matching point coordinate 


x* + x * 


A^mi, 


(A5) 


and then from (74) 


a > a i ~ 2 ( X ~ e 2X * 1 ) ■ 


(A6) 


Equations (A2), (A4), ( |A5[ ) and ( |A6[ ) define a\ as a func¬ 
tion of A x , such that for a > ai(A^) there can be more 
than one discrete solution corresponding to the same a. 
The graph of this function is shown in fig. |19| Similarly, 
by solving the inequality 


®m +2 


< a n 


( n 3 + ji~ 3 

\ (M + 1) ’ 

h \ ^2m +1 _ ^ j 

{ (fi + 1 ) 


j £ 0,TO+ 1, 

(A3) 

j £ TO, OO. 


This result is valid under the assumption that 


& £ (Um+1 > Urn) 

= (a m , a . m ). 


///" - fj,-™- 1 / u m+1 - \ 

V M m+1 + ’ M m+1 + J 


we get the function a 2 {A x ) such that for a > a 2 (A X ) 
the non-uniqueness of the discrete solutions is not only 
possible, but is guaranteed (we omit the straightforward 
but bulky derivation). The graph of this function is also 
shown in fig. [To] 

However the question of stability of the discrete critical 
nucleus solution is more important, even if this solution 
is unique. For Uj(t) = Uj +Vj(t), the linearized system is 


dv j 
d t 


Vj- 1 - 2 Vj + Vj+ 1 


u 3> 


j £ 


(A7) 


This gives a range of possible values of a for a given 
to. So, at a fixed A x , the dependence of the solution on 
parameter a is discontinuous (piecewise constant), and 
there is a possibility that at some combinations of A x 
and a, there could be more than one solution. Indeed, 
this possibility is realized if the intervals ( a m ,a m ) for 
consecutive values of to overlap, that is, 


®m+l 

which is the case whenever 

log (yi 2 + n + l) 3 


TO > TOi = 


2 log /T 


(A4) 


when one of Uj = a exactly. The spectrum of the sys¬ 
tem (A7) in £ 2 (Z) is [—1 — 4/A x 2 , —1], with eigenpairs 


Wj = exp (inj), X = — 1 — 2(1 — cosk)/A 




K £ 


So the discretized critical nuclei are almost surely asymp¬ 
totically stable in the linear approximation (remember 
that we did not analyse the exceptional cases when 
Uj = a exactly at some j). As a result, in our sim¬ 
ulations at initiation parameters close to the threshold, 
we find distinct “frozen critical nucleus” solutions, so that 
the critical curve between ignition and failure becomes in 
fact a “critical band”. “Convergence” of the discretized 
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FIG. 20: (Color online) The region where “frozen nuclei” solutions occur in direct numerical simulations of the McKean 
model ( |72| ) for a = 0.45 and x s = 0.5 and with variation of the space step A x is located between the two dotted red lines with 
annotated data points and is shaded in gray. The thin solid blue lines are lines of linear extrapolation the intersection of which 
indicated by the violet star and annotated is taken as the boundary between decay and ignition, (b) Comparison between the 
“frozen” solution obtained by direct numerical simulation at a = 0.45, x s = 0.5, U s = 2.13, A x = 0.01 and the known exact 
analytical expression for the critical nucleus (|73|) at o = 0.45. 


system to the continuous system, as far as the initia¬ 
tion problem is concerned, is manifested by reduction of 
the basin of attraction of the critical nucleus solutions 
with decreasing discretization steps, as is also evident in 
fig. [19] This convergence is illustrated in fig. [20| 


Appendix B: On linearization of discontinuous 
right-hand sides 


Rinzel and Keller [42] obtained operator (75) by for¬ 
mally differentiating the nonlinearity of (72), apparenly 
having in mind a calculation like 


Let us consider a one-parametric family of solutions, 
u = f/(x,t;p), which is continuous and piecewise- 
differentiable in a;, t and p. We are seeking a differential 
equation for v = dU/dp , that is, the linearized equation, 
or equation in variations. Since the function v is expected 
to be discontinuous, its differential equation will contain 
singular terms, i.e. it should be understood in terms of 
generalized functions (distributions). 

Suppose that our solution is monotonically decreasing 
in x in some region of the ( x , t) plane containing the curve 
x = s(t) at which U changes sign (the opposite case is 
considered in a similar way). Then (B1) can be rewritten 
as 


df(u) df(u(x)) / du(x) 


du 


dx 

-l 


dx 


fdu(x)\ d , , 


1 


u'(x) 


( — u'{x) - 5(x* — x)) 
S(x - x*). 


•i'(x*) 


This might look paradoxical, as the linearization pro¬ 
cedure in its traditional understanding is based on the 
assumption of smallness of the increments, whereas dis¬ 
continuity of the reaction term means that in some cir¬ 
cumstances increments of those terms are not small when 
the increments of the arguments vanish, which is an ap¬ 
parent contradiction. Nevertheless, this procedure can 
be mathematically justified. 

Consider, as a simple example, a scalar reaction- 
diffusion equation for some field u(x,t ), with a single 
discontinuity at u = 0: 


(Bl) 


Ut = h(U) + (f 2 (U) - A([/))H (s(t) -x) + DU XX , 
U(s(t),t;p) = 0. 

By differentiating this with respect to p, we have 

Ut P = fi'(U)U p + ( h'(U ) - h\u)) H(s(i) - X )U P 
+ DU XX p + (M0) - M0)) 8(s(t) - x) s p , 
U p (s(t),t;p ) + U x (s(t),t;p) s p = 0. 

By excluding s p from this system and setting U p = v, we 
have 


Vt = fi(U)v + (M(U) - fi(U)) H(s(f) - x)v 

M0)-M0) 


Dv x 


U x (s,t\p) 


■ S(s(t) — x) v, 


which is precisely the linearized equation which would be 
obtained by formal differentiation of the right-hand side 
of (Bl), and subsequent replacement of S(U) by 5(s — x) 
using the chain rule. 


u t = M u ) + (Mu) - fi(u))R(u) + Du, 
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Appendix C: Finite element discretization for the 
McKean model 


For the McKean model, we needed to solve initial-value 
problems both for the nonlinear equation, 


du d 2 u 
dt = dx? 

and its linearization, 

dv d 2 v 


— u + H(it — a) 


— v + S(u — a)v. 


dt dx 2 

These are equivalent (for u decreasing in x ) to 
du d 2 u 

m = + 

dv d 2 v 1 . 

= 7W ~ v + “77 —\ d \ x - x *) v 
dt ox 1 u'[x *) 

where x* = x»(f) is defined by 

n(x*(f), t) = a. 

Both cases required finite-element treatment, and we 
present the details for both cases together, by writing 
both as 

dw d 2 w 

Si = SS + IM 

where w = u, f = —u + H(it — a) = —it -(- H(x* — x) 
in one case, and w = v, f = — v + S(u — a)v = —u + 
(l/u'(x*))8(x — x*)i> in the other case. 

The finite element method (see e.g. [55]) is based on a 
weak formulation of the problem, which requires that 

l 4>w {^“5^“ /( ”’’ iE) } da: = 0 

for any “test function” 4>(x). If the variety of the avail¬ 
able test functions is broad enough then the weak for¬ 
mulation is equivalent to the original equation. After 
integration by parts and taking into account the Neu¬ 
mann boundary conditions for w, the weak formulation 
can be formally re-written as 


, (dw 

* w (a - 


ff u , 90 dw 


dx = 0. (Cl) 


The difference is, of course, that whereas the original 
formulation requires that w is twice differentiable in x, 


the weak formulation (Cl) uses only first derivatives of 


w, and can be applied as long as the test functions $ are 
once differentiable. 

The standard finite element method is the Galerkin 


method applied to (Cl). It uses a set of linearly inde¬ 
pendent functions, $J[x) , j = 0,..., N, called the finite 


elements, as the test functions, and seeks the approxima¬ 
tion of the solution in the span of this same set: 


N 


'j(x, t) ~ 1/7(x, t) = Wj(t)<&j(x). (C2) 

7=0 


Substitution of (C2) into (Cl) for $ = i = 0,..., N, 


leads to the system of equations 
N ^ N 

^ ' -Aij ^ ^ ' BijWj = Ci ( uij ), * = 0,..., N, 

1=0 1=0 

or in the vector form, for w (t) = 


A A +B w = C 

dt 

where the coefficients are given by 


A%,j — 


4>,:(x)<F,(x) dx, 


Bi,j — 


) dx, 


(C3) 

(C4a) 

(C4b) 


f 4>'(x)4>'(x) 

Jo 

Ci(wj)= J ®i{x)f dx. (C4c) 

We use a simple and popular choice of the test functions, 
the piecewise linear tent functions: 

{ (x - x<_i) /A x , xG[xi_i,Xi], 

(xj+i - x) /A x , x G [xi, x i+ i], (C5) 

0, otherwise, 

for a regular grid of (xj), 

Xi =i A x , i = 0,..., N, A X =L/N. (C6) 

Obviously, in this case <F,(xi) = Sj , 7 - and therefore 
w(xi) = vbi. For these test functions, (C4) give the mass 

matrix A = 


as 


A x 
A = — 
6 


(2 1 0 0 \ 

14 1 : 

o 0 

: 14 1 

Vo ••• 0 12/ 


(C7) 


as 


the stiffness matrix B = 

( 1 -1 0 
-1 2 -1 
0 


1 


B = — 


0 \ 


0 


V 0 


-1 2 -1 
0 -1 1 / 


(C8) 
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as 


and the load vector C = 

C (w) = —Aw + F, 


(C9) 


where F = (j?ij and differs for the nonlinear problem 
and for the linearized problem. 

For the nonlinear problem, w = u, we get F = F^ 1 ) + 
F^ 2 ), where 


Ui -1 > a, v,i > a, 
Ui-i > a, Ui < a, 
Ui -1 < a, Ui> a, 
otherwise, 

(CIO) 


^2+1 V*i 

Ui+1 > a, Ui < a, 

5 ^2+1 ^ ^2 

otherwise, 

(Cll) 



and 


F (2) = 1 


2A, 


A,, 2 , 

(x* - x i+ i) 2 , 

A x (x* x i+ i) 

. 0 , 


for z = 1,..., A — 1, and 
f2A, 2 , 


Fn = 


1 


4A, 


and 


F n = 


4A a 


(x* - X_i) 2 
+ (xi - x») 2 , 

2A X 2 - (x* - x_i) 2 
- (xi - x*) 2 , 

lo, 


f2A* 2 , 

(x* - Xat+i) 2 
+ (xjv-i - x*) 2 , 

2A X (x* xat-j-i) 
- (xjv-1 - X*) 2 , 


v0, 


uq > a, u i > a, 

ilq < a, u\ > a, 

Uq > a, Ui < a, 
otherwise, 


UjV—1 > CL, Um > CL, 

un -i > n, ww < a, 

UN-1 < d, dw > Cl, 
otherwise. 

(C13) 


for the boundary points. In these formulas, x» is the 
point such that w(x*,f) = a by linear interpolation, i.e. 
for m such that (u m — a)(u m +1 — a) < 0, we have x* = 
((«m+i - a )Xm + (d - d m )x m+ i) /(u m+ i - u m ), and the 
definition (C6| is extended to i = — 1 and i = N + 1. 


For the linear problem, w = v, we get 

1 


p — 

1 m. — 


X*) Uri 


aA x 2 

T (*^m+1 *£*) (*£* ■Z'm) fim+l] 


F m +1 — 


1 

aA x 2 


(C14) 

[(Xm+l X*) (x* x m ) 

T (x* x m ) x m +i , (C15) 

Fj =0, j ^ m,m+1, (C16) 

where m and x* are defined based on the nonlinear solu¬ 
tion u based on the same rule as above. 

Appendix D: The modified Beeler-Reuter model 


The model was proposed in 


We use it in 

,T 


the following formulation: f : (V,h,j,X\,d,f,Ca) i-» 

ifv, fh , fj 5 f xi , fd, ff, /ca) T , where 

fv = —I Kl fY) ~ Ixi (V) xr) 

-I Na (V,h,j) - I s (V,d, f, Ca), 
h = a h (V)(l-h)-p h (V)h, 
fj = a 0 {V)(l-j)-pj(V)j, 

fx 1 = CX X i(F)(l - Xi) - ^i(F)xi, 

fd = a d (V)(l - d) - /3d(V)d, 
f f = a f (V)(l-f)-p f (V)f, 

/Ca = —10 _7 / s + 0.07(10~ 7 — Ca), 


(C12) where the transmembrane currents are defined by 


Ik i {V) = 0.35 (0.3 — a) I Kl (V), 

_ 4 t e 0.04(V+85) _ i') 

IkAV)= { ’ 


e 0.08(y+53) _|_ e 0.04(y+53) 

0.2(V + 23) 


1 _ e -0.04(y+23) ’ 

I x l(V, Xl ) — C/ix( f 7 )xi, 

e 0.04(y+77) _ I 
9ix{V) — 0.8 e 0.04(y+35) ’ 

lNa(V,h,j) = (g Na ( m(V)f hj + g Nac )(V - E Na ), 
Is(V,d, f, Ca) = g s df{V - E s { Ca)), 

the rn-gate is assumed in quasi-stationary state, 

m(V) = a m (V)/(a m (V) + j3 m (V)), 
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the gate opening and closing rates are defined by 
0.0005 e°- 083 ( y + 50 ) 


a xl (V) = 


Pxi(V) = — 
Om(V) = 


e 0.057(y+50) I l ’ 

0.0013 e-° 06 ( y+20) 


e -0.04(V+20) _|_ l 

V + 47 

X _ e -0.1(V+47) ’ 


/3 m (V) =40e-° 056(y+72) , 
a h (V) = 0.126e“°- 25(y+77) , 


i3 h (V) = - 


1.7 


e -0.082(^+22.5) _|_ X ’ 


i{V) = — 


0.055 e -°- 25 ^+ 78 ) 


0.2(y+78) _)_ X 


Pj(V) 

a d (V) 


0.3 

e -0.1(V+32) + x ’ 

0.095 e -0-Oi(v-5) 

g-0.072(V-5) I X ’ 


MV) 

a f (V) 

Pf(V) 


0.07e _ool7 ( y+44 ) 

g0.05(V+44) _|_ X ’ 

0.012 e -o.oo 8 (v+ 28 ) 

e 0.15(V+28) _|_ X 

0.0065 e-°- 02 ( y + 3 °) 


e -0.2(V+30) _)_ X 


and the calcium reversal potential is defined by the 
Nernst law, 


E s { Ca) = -82.3 - 13.0287log(Ca). 


The parameters of the model are fixed at the values used 
in I47H50] : g Na = 0.24, g Nac = 0.003, E Na = 50, g s = 
0.045, and for a we consider two values a = 0.105 and 
a = 0.115. 
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